{"cells":[{"cell_type":"markdown","metadata":{"editable":true,"id":"A51DDC3D042E44A5B2C11E681044CFB2","jupyter":{},"mdEditEnable":false,"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":["figure"]},"source":[" # <center> Lecture 3: The Beta-Binomial Bayesian Model </center>  \n"," \n"," ## <center> Instructor: Dr. Hu Chuan-Peng </center>  \n"]},{"cell_type":"markdown","metadata":{"id":"9B220185A8284B738BFF11F846D9D12A","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":[" ## <center> Part 1: Beta 先验 </center>"]},{"cell_type":"markdown","metadata":{"editable":true,"id":"228798C5F4F24D329CF8BC1A7D3096FF","jupyter":{},"mdEditEnable":false,"notebookId":"66ea43e4f99628c505715aea","notebookRunGroups":{"groupValue":"1"},"runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":["figure"]},"source":["上节课的内容主要是关注了贝叶斯的思路：即我们如何根据信息更新自己对某个事件/变量的信念，并顺便比较了一下贝叶斯和频率学派对于同一事件的不同解读视角。  \n","\n","🎯接下来，我们将贝叶斯理念融入心理学研究中，探索其如何优化我们对认知过程的推理过程。"]},{"cell_type":"markdown","id":"b1b5a47b","metadata":{"editable":true,"id":"A63998DCB1914A7DB9870477E5E1FFD1","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["### 以随机点运动任务（Random motion dot task）为例  \n","\n","- 随机点运动任务是一个经典的认知任务，是研究感知决策(perceptual decision making)的标准范式。  \n","- 被试任务：  \n","   - 被试观察屏幕上的点，这些点以不同的一致性水平随机移动。  \n","   - 任务是判断这些点的主要移动方向（例如，向左还是向右）。  \n","- 实验可以在人类和动物（如老鼠）上进行，以探索人类和动物的大脑如何完成这一过程。  \n","\n","<center>  \n","    <table>  \n","            <tr>  \n","                <td><img src=\"https://cdn.kesci.com/upload/sjwnyi477j.gif?imageView2/0/w/400/h/400\" alt=\"\"></td>  \n","                <td><img src=\"https://cdn.kesci.com/upload/sjwnyt1yq4.gif?imageView2/0/w/400/h/400\" alt=\"\"></td>  \n","            </tr>  \n","            <tr>  \n","                <td>一致性10%</td>  \n","                <td>一致性60%</td>  \n","            </tr>  \n","    </table>  \n","</center>  \n","\n","-----------  \n","-----------"]},{"cell_type":"markdown","id":"58987e6a","metadata":{"editable":true,"id":"EFDA064DF3B0418992FCD1C7360B0EE9","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["正确率与任务难度：  \n","   - 任务的难度由点的运动一致性决定。一致性低，点的主要方向不明显，任务难度大；一致性高，点的主要方向明显，任务难度小。  \n","   - 正确率随一致性的增加而提高，通常在心理物理曲线中展示。  \n","\n","根据正确率的数据，我们可以使用心理物理曲线描述运动强度（点的一致性百分比）与决策正确率之间的关系。  \n","\n","![image](https://cdn.kesci.com/upload/sjzjcgwjrm.png?imageView2/0/w/640/h/640)  \n","\n","> Shooshtari, S. V., Sadrabadi, J. E., Azizi, Z., & Ebrahimpour, R. (2019). Confidence representation of perceptual decision by EEG and eye data in a random dot motion task. Neuroscience, 406, 510–527. https://doi.org/10.1016/j.neuroscience.2019.03.031  \n","\n","------  \n","\n","------"]},{"cell_type":"markdown","id":"ee9284bf","metadata":{"editable":true,"id":"6E0BE9840A634ADA9E06FE9C95643BEF","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["\n","正确率这个指标与被试在每个试次上的关系是什么？回想一下刚刚讲到的频率主义与贝叶斯主义的区别。  \n","\n","正确率给出一个单一的指标，而不是对无法观察的“被试做正确的能力”进行估计，也没有纳入估计的不确定性。  \n","\n","假如从贝叶斯的角度出发，我们根据被试在某种条件下的表现来估计其做正确的能力？  \n","\n","我们以Evans et al.（2020, Exp. 1） 的数据为例进行探索。  \n","\n","> Evans, N. J., Hawkins, G. E., & Brown, S. D. (2020). The role of passing time in decision-making. Journal of Experimental Psychology: Learning, Memory, and Cognition, 46(2), 316–326. https://doi.org/10.1037/xlm0000725  \n"]},{"cell_type":"code","execution_count":1,"id":"2178953a","metadata":{"collapsed":false,"id":"512CD3D6D1E34D8E9A10ACC6C65DC7F8","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","scrolled":false,"slideshow":{"slide_type":"skip"},"tags":[],"trusted":true},"outputs":[],"source":["# 导入必要的库\n","import scipy.stats as st\n","import numpy as np\n","import pandas as pd\n","import seaborn as sns\n","import matplotlib.pyplot as plt\n","import preliz as pz\n","\n","from scipy.stats import binom\n","\n","# 为 preliz 绘图设置图形样式\n","pz.style.library[\"preliz-doc\"][\"figure.dpi\"] = 100\n","pz.style.library[\"preliz-doc\"][\"figure.figsize\"] = (10, 4)\n","pz.style.use(\"preliz-doc\")\n","\n","# 使用 pandas 导入示例数据\n","try:\n","  data = pd.read_csv(\"/home/mw/input/bayes3797/evans2020JExpPsycholLearn_exp1_clean_data.csv\") \n","except:\n","  data = pd.read_csv('data/evans2020JExpPsycholLearn_exp1_clean_data.csv')\n","\n","# 选择特定模式的数据进行演示：\n","# 筛选 correct 为 1 的数据，并随机选取30个\n","data_correct_1 = data[data['correct'] == 1].sample(n=30, random_state=42)\n","\n","# 筛选 correct 为 0 的数据，并随机选取20个\n","data_correct_0 = data[data['correct'] == 0].sample(n=20, random_state=42)\n","\n","# 合并两部分数据，形成新的数据集\n","df = pd.concat([data_correct_1, data_correct_0], ignore_index=True)"]},{"cell_type":"code","execution_count":2,"id":"fe0fc7ab","metadata":{"collapsed":false,"id":"FBF5B2F9E81B498DAAF90C49860AEFDC","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>subject</th>\n","      <th>numberofDots</th>\n","      <th>percentCoherence</th>\n","      <th>correct</th>\n","      <th>RT</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>0</th>\n","      <td>82111</td>\n","      <td>40</td>\n","      <td>5</td>\n","      <td>1</td>\n","      <td>1173</td>\n","    </tr>\n","    <tr>\n","      <th>1</th>\n","      <td>82111</td>\n","      <td>40</td>\n","      <td>5</td>\n","      <td>1</td>\n","      <td>1820</td>\n","    </tr>\n","    <tr>\n","      <th>2</th>\n","      <td>82111</td>\n","      <td>40</td>\n","      <td>5</td>\n","      <td>1</td>\n","      <td>694</td>\n","    </tr>\n","    <tr>\n","      <th>3</th>\n","      <td>82111</td>\n","      <td>40</td>\n","      <td>5</td>\n","      <td>1</td>\n","      <td>1355</td>\n","    </tr>\n","    <tr>\n","      <th>4</th>\n","      <td>82111</td>\n","      <td>40</td>\n","      <td>5</td>\n","      <td>1</td>\n","      <td>1480</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["   subject  numberofDots  percentCoherence  correct    RT\n","0    82111            40                 5        1  1173\n","1    82111            40                 5        1  1820\n","2    82111            40                 5        1   694\n","3    82111            40                 5        1  1355\n","4    82111            40                 5        1  1480"]},"execution_count":2,"metadata":{},"output_type":"execute_result"}],"source":["# 只显示特定列\n","df[['subject', 'numberofDots', 'percentCoherence', 'correct', 'RT']].head()"]},{"cell_type":"code","execution_count":3,"id":"dc306b1b","metadata":{"collapsed":false,"id":"B2FAFFF5E0534577B09443174E8A620E","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"outputs":[],"source":["# 非选课的同学，可以使用和鲸社区的“2024社区镜像”，运行以下的代码安装必要的模块，需要3-5分钟的时间完成加载\n","# 后续会有专门的社区公开镜像，给大家提前配置好运行环境\n","# 将下列代码行解除注释，删除“#”，运行即可：\n","# !conda install -y graphviz bambi=0.13.0 pymc=5.16.2 PreliZ=0.9.0 ipympl=0.9.4 pingouin=0.5.4\n","\n","# docker部署和使用教程链接：https://zhuanlan.zhihu.com/p/719739087\n","# docker pull hcp4715/pybaysian:latest\n","# docker run -it --rm -p 8888:8888 hcp4715/pybaysian:latest"]},{"cell_type":"markdown","id":"f1626b47","metadata":{"id":"1899A56BDFC14D25B8A9D4FD29FB8FC0","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["##### 接下来，我们将通过该数据学习以下知识点：  \n","\n","1. 如何为结果为“正确与否”的数据建立概率模型  \n","2. 如何为该类数据设置一个先验模型  \n","3. 如何通过数据来更新先验模型"]},{"cell_type":"markdown","id":"65f657d9","metadata":{"editable":true,"id":"A3E6CEA11DB94A32AD21352B93625B8E","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["### 连续变量的先验模型"]},{"cell_type":"markdown","id":"181f70e6","metadata":{"editable":true,"id":"0145FB1B0CFA43C1851CDECA802EB7F5","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":""},"tags":[]},"source":["#### \"正确与否(0/1)数据\"的处理  \n","\n","在随机点运动任务中，被试正确探测出点的运动方向的取值范围在0-1之间，不同被试、不同条件下的真实的能力不同（与上节课的“重复成功率”类似）。  \n","\n","假如我们对人们在5%一致条件下进行正确判断的能力感兴趣。  \n","\n","- 如果我们还没有招募被试进行实验，则只能对被试在不同条件下的表现进行推断。根据先前的文献(Vafaei Shooshtari et al., 2019)，当随机点的一致性为5%时，个体的正确率约为70%  \n","\n","  > Shooshtari, S. V., Sadrabadi, J. E., Azizi, Z., & Ebrahimpour, R. (2019). Confidence representation of perceptual decision by EEG and eye data in a random dot motion task. Neuroscience, 406, 510–527. https://doi.org/10.1016/j.neuroscience.2019.03.031  \n","\n","- 根据文献的信息，我们可以建立一个关于正确判断能力的先验分布：例如，个体正确率为70%的概率为0.8，而个体正确率为60%或80%的概率为0.1 (这是我们上节课学习的内容🤓)。  \n","\n","但考虑到正确率为0-1之间的任意值，我们只采用0.6、0.7和0.8三个值，似乎有点武断。可能连续概率模型(continuous prior probability)更合适?  \n","\n","🤔如何建立连续变量的先验概率模型？  "]},{"cell_type":"code","execution_count":4,"id":"a7f46e8b","metadata":{"collapsed":false,"id":"BC3FFDAD60D6433D8F16F9989140CC9D","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/BC3FFDAD60D6433D8F16F9989140CC9D/sjznu3v1zc.png\">"],"text/plain":["<Figure size 800x500 with 1 Axes>"]},"metadata":{},"output_type":"display_data"}],"source":["# 创建离散先验分布数据\n","prior_disc = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]         # 离散先验分布的取值\n","prior_dis_prob = [0, 0, 0, 0, 0, 0.10, 0.80, 0.10, 0, 0]  # 对应的概率\n","prior_disc_data = pd.DataFrame({'ACC': prior_disc, 'f(ACC)': prior_dis_prob})\n","\n","# 将'ACC'列的数据类型转换为字符串，以便在图表中正确显示\n","prior_disc_data['ACC'] = prior_disc_data['ACC'].astype('str')\n","\n","# 创建一个单独的图表\n","plt.figure(figsize=(8, 5))\n","\n","# 绘制离散先验分布的条形图\n","sns.barplot(data=prior_disc_data, x='ACC', y='f(ACC)', palette=\"deep\")\n","\n","# 设置图表的标题\n","plt.title(\"Discrete Prior\")\n","\n","# 移除图的上边框和右边框\n","sns.despine()\n","\n","# 显示图表\n","plt.show()"]},{"cell_type":"markdown","id":"1dfcfe16","metadata":{"id":"3BA73CF8568D4ED6BD242899129BB041","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["### Beta 先验模型"]},{"cell_type":"markdown","id":"d5389b93","metadata":{"id":"3463E2AAD2414A659E42CCF87BD534D4","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["⏰回顾一下离散型随机变量及其**概率质量函数(probability mass function, pmf)**  \n","\n","> 在随机点运动任务的例子中，假如正确率(ACC)只有几个可能的取值，则其是离散型随机变量$Y$，可以通过概率质量函数来描述离散型随机变量在各特定取值上的概率，对所有正确率的取值来说，$0\\leq ACC \\leq 1$， 并且 $\\sum_{all\\,\\pmb{y}}f(ACC) = 1$，ACC取值的所有概率之和为1。"]},{"cell_type":"markdown","id":"9cc122ba","metadata":{"id":"90846BA613AF448A914126B9FCE0D2A7","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["**Continuous probability models**  \n","\n","当正确率是一个连续的随机变量时，它服从于一个分布$f(ACC)$, 又被称概率密度函数(probability density function, pdf)  \n","\n","pdf 与 pmf 有类似的性质：  \n","\n","- $f(ACC) \\ge 0$  \n","- $\\int_{\\text{ACC}} f(\\text{ACC}) \\, d\\text{ACC} = 1$，即 $f(\\text{ACC})$ 曲线下的面积之和为 1  \n","- 当$a \\le b$时，$P(a < ACC < b) = \\int_a^b f(ACC) dACC$  \n","\n","![Image Name](https://cdn.kesci.com/upload/s0ywgsar31.jpg?imageView2/0/w/500/h/500)  \n","> source:https://www.zhihu.com/question/263467674/answer/1117758894  \n","\n","--------------  \n","--------------"]},{"cell_type":"markdown","id":"040104f2","metadata":{"id":"47A218285789430F90A7BA1B7412D9A7","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["**为正确率（ACC）选择一个合适的分布 -- Beta分布**  \n","\n","虽然我们潜意识里可能认为正确率服从正态分布，但是，在统计世界中存在着许多分布。  \n","\n","在这里，我们假设随机点运动任务中的正确率服从Beta分布。  \n","\n","- 选择 Beta 分布的关键在于，其取值范围满足[0,1]，  \n","- 此外，Beta分布通过两个超参数(hyperparameters)，$\\alpha \\; (\\alpha>0)\\;和 \\beta\\;(\\beta>0)$ 调节分布的形态。  \n","\n","$$  \n","\\text{ACC} \\sim \\text{Beta}(\\alpha, \\beta)  \n","$$  \n"]},{"cell_type":"markdown","id":"5c795972","metadata":{"id":"B5761E61DB55409D97646C51F9B44761","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["**Beta分布的参数：控制分布形状—— $\\alpha$和$\\beta$**  \n","\n","调整$\\alpha$和$\\beta$的大小，观察分布形状的变化  \n","\n","下图展示了不同$\\alpha$和$\\beta$下pdf形态的变化：  \n"]},{"cell_type":"code","execution_count":5,"id":"088df461","metadata":{"collapsed":false,"id":"E2B69AFCDB0A47B28F33857C985213DC","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/E2B69AFCDB0A47B28F33857C985213DC/sjznu79qv7.png\">"],"text/plain":["<Figure size 1000x1000 with 9 Axes>"]},"metadata":{},"output_type":"display_data"}],"source":["# 创建一个3x3的网格子图，每个子图的尺寸为10x10\n","fig, axs = plt.subplots(3, 3, figsize=(10, 10))\n","\n","# 绘制 Beta(1, 5) 分布的PDF，并显示置信区间\n","pz.Beta(1, 5).plot_pdf(pointinterval=True,ax=axs[0, 0], legend=\"title\")\n","pz.Beta(1, 2).plot_pdf(pointinterval=True,ax=axs[0, 1], legend=\"title\")\n","pz.Beta(3, 7).plot_pdf(pointinterval=True,ax=axs[0, 2], legend=\"title\")\n","pz.Beta(1, 1).plot_pdf(pointinterval=True,ax=axs[1, 0], legend=\"title\")\n","pz.Beta(5, 5).plot_pdf(pointinterval=True,ax=axs[1, 1], legend=\"title\")\n","pz.Beta(20, 20).plot_pdf(pointinterval=True,ax=axs[1, 2], legend=\"title\")\n","pz.Beta(7, 3).plot_pdf(pointinterval=True,ax=axs[2, 0], legend=\"title\")\n","pz.Beta(2, 1).plot_pdf(pointinterval=True,ax=axs[2, 1], legend=\"title\")\n","pz.Beta(5, 1).plot_pdf(pointinterval=True,ax=axs[2, 2], legend=\"title\")\n","\n","# 自动调整子图之间的间距，以防止标签重叠\n","plt.tight_layout()\n","\n","# 显示绘制的图形\n","plt.show()"]},{"cell_type":"markdown","id":"ac08c534","metadata":{"id":"D11DE09CF8E9408B9B2C583219B98DFB","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["**Beta分布的概率密度函数pdf**  \n","\n","$$  \n","f(\\text{ACC}) = \\frac{\\Gamma(\\alpha + \\beta)}{\\Gamma(\\alpha)\\Gamma(\\beta)} \\text{ACC}^{\\alpha-1} (1-\\text{ACC})^{\\beta-1} \\;\\;  \\text{for} \\; \\text{ACC} \\in [0,1]  \n","$$  \n","\n","*简单了解：*  \n","\n","* *$\\Gamma$与阶乘有关*  \n","\n","* *$\\Gamma(z) = \\int_0^\\infty x^{z-1}e^{-y}dx$ 且 $\\Gamma(z + 1) = z \\Gamma(z)$*  \n","\n","* *当z是正整数的时候，$\\Gamma(z)$ 可以被简化为$\\Gamma(z) = (z-1)!$*"]},{"cell_type":"markdown","id":"7437529c","metadata":{"id":"2BEA2C2C5B68424991C60264045D5FFD","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["**调整Beta先验**  \n","\n","在这里我们选择$\\text{ACC} \\sim \\text{Beta}(70, 30)$作为合理的先验模型  \n","\n","* *关于为何$\\alpha=70, \\beta=30$，感兴趣的同学可自行查看课后bonus💐*  \n","\n","带入公式，可以计算出先验$f(ACC)$的pdf：  \n","\n","- pdf:  \n","$$  \n","f(ACC) = \\frac{\\Gamma(100)}{\\Gamma(70)\\Gamma(30)} ACC^{69}(1-ACC)^{29} \\;\\; \\text{ for } \\; ACC \\in [0,1] .  \n","$$  \n"]},{"cell_type":"markdown","id":"3fb52850","metadata":{"id":"5C7AC8C1DF184B2B9C11E42AEAF40A47","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["以下就是我们这个先验似然的可视化： "]},{"cell_type":"code","execution_count":6,"id":"1eec8a46","metadata":{"collapsed":false,"id":"C6445573C12A44019380BA71B4E4908C","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/C6445573C12A44019380BA71B4E4908C/sjznu73qas.png\">"],"text/plain":["<Figure size 1500x400 with 2 Axes>"]},"metadata":{},"output_type":"display_data"}],"source":["# 创建离散先验分布数据\n","prior_disc = [0.1, 0.2, 0.3,0.4,0.5,0.6, 0.7, 0.8,0.9, 1.0]         # 离散先验分布的取值\n","prior_dis_prob = [0, 0, 0, 0, 0, 0.10, 0.80, 0.10, 0, 0]  # 对应的概率\n","prior_disc_data = pd.DataFrame({'ACC': prior_disc, 'f(ACC)': prior_dis_prob}) \n","# 将'ACC'列的数据类型转换为字符串，以便在图表中正确显示\n","prior_disc_data['ACC'] = prior_disc_data['ACC'].astype('str')\n","\n","# 创建一个包含两个子图的图表\n","f, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 4))\n","\n","# 在第一个子图中绘制离散先验分布的条形图\n","sns.barplot(data=prior_disc_data, x='ACC', y='f(ACC)', palette=\"deep\", ax=ax1)\n","# 设置第一个子图的标题\n","ax1.set_title(\"discrete prior\")\n","\n","# 在第二个子图中绘制连续先验分布的线图\n","pz.Beta(70, 30).plot_pdf(pointinterval=True,ax=ax2, legend=\"title\") # 在这里示范了一个 Beta 分布的参数，你可以根据需要修改这些参数\n","# # 设置第二个子图的标题\n","ax2.set_title(\"continuous prior\\n(Beta(alpha=70,beta=30))\")\n","\n","# 设置第二个子图的 x 轴范围为 0 到 1\n","ax2.set_xlim(0, 1)\n","\n","# # # 移除图的上边框和右边框\n","sns.despine()"]},{"cell_type":"markdown","id":"9cc373b7","metadata":{"id":"AC842485216B4C01993B90B8BBDCE312","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["## <center>Part 2 Binomial 与似然 </center>"]},{"cell_type":"markdown","id":"0c0e8a01","metadata":{"id":"FE24090F27B34510954216BCE485DAF6","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["### 二项数据模型&似然函数"]},{"cell_type":"markdown","id":"7a8a51ad","metadata":{"id":"44D327A99B514238A6DCF328C640DC14","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["\n","确定先验之后，第二步，我们需要根据收集的数据构建似然函数  \n","- 我们选取了 Evans 等(2020)的研究数据作为我们的案例展示，为计算方便，我们仅选取一名被试的数据，并对该被试数据进行选择以方便展示。  \n","- 假设该被试的数据包含 50 次试验的随机点运动任务，我们希望通过 50 次的实验来估计被试在5%一致性条件下的判断正确的能力  \n","  - 我们假设参与者的每一次判断是否正确是相互独立的。  \n","- 假设参与者的正确判断次数为 Y，我们可以做出以下假设：  \n","\t1. 每次试验的正确与否是相互独立的（试验A的结果不会影响试验B）。  \n","\t2. 每次试验参与者正确判断的概率为ACC。"]},{"cell_type":"markdown","id":"05f0bbc2","metadata":{"id":"1F96853F478446B7AE54F534616459CD","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["在这两个假设下，正确判断次数 Y 和正确率ACC之间的关系符合二项分布：  \n","$$  \n","Y | \\text{ACC} \\sim \\text{Bin}(50, \\text{ACC})  \n","$$  \n","\n","在不同的正确率下，出现特定正确判断次数Y的概率$f(y|\\text{ACC}) \\quad y \\in \\{0,1,2,…,50\\}$ 可以表示为：  \n","\n","$$  \n","f(y|\\text{ACC}) = P(Y=y | \\text{ACC}) = \\binom{50}{y} \\text{ACC}^y (1-\\text{ACC})^{50-y}  \n","$$"]},{"cell_type":"markdown","id":"e7950b6f","metadata":{"id":"AE3336614ECE4F27BC738FD28FC934ED","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["我们来进行一些具体的计算  \n","- 假设ACC = 0.1，那么在50 次试验中，每个可能的正确判断次数$Y \\in ({1,2...50})$，而出现特定结果的概率可以表示为：  \n","$$  \n","f(Y=1|\\pi=0.1) = \\binom{50}{1} 0.1^1 (1-0.1)^{49}  \n","$$  \n","$$  \n","f(Y=2|\\pi=0.1) = \\binom{50}{2} 0.1^2 (1-0.1)^{48}  \n","$$  \n","$$  \n","...  \n","$$  \n","$$  \n","f(Y=49|\\pi=0.1) = \\binom{50}{49} 0.1^{49} (1-0.1)^{1}  \n","$$  \n","$$  \n","f(Y=50|\\pi=0.1) = \\binom{50}{50} 0.1^{50} (1-0.1)^{0}  \n","$$  \n","\n","我们可以把这50个概率值画出来。"]},{"cell_type":"code","execution_count":7,"id":"aabd3e25","metadata":{"collapsed":false,"id":"E380B66D21C44094987B4A66A4E1FED5","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/E380B66D21C44094987B4A66A4E1FED5/sjznu839xh.png\">"],"text/plain":["<Figure size 1000x400 with 1 Axes>"]},"metadata":{},"output_type":"display_data"}],"source":["# 创建一个 Binomial 分布，参数为 n=50（试验次数），p=0.1（正确率概率）\n","binom = pz.Binomial(n=50, p=0.1)\n","\n","# 绘制该分布的概率密度函数 (PDF)\n","binom.plot_pdf()\n","\n","# 移除图的上边框和右边框\n","sns.despine()\n","\n","# 显示图表\n","plt.show()"]},{"cell_type":"markdown","id":"993d3af3","metadata":{"id":"0D8E72C31217475A9058931EB6CBC959","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["**随机点运动任务中的正确率分布**  \n","\n","在随机点运动任务中，正确率描述了参与者正确判断点运动方向的概率。在不同的正确率下，参与者在50次实验中正确判断的次数Y会有不同的分布情况。  \n","\n","我们可以模拟不同的ACC，并展示每个ACC下，正确判断次数Y的分布图。  \n","\n","（注：这里只展示了当正确率 =0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 九种情况下的分布图，但是记住ACC的取值其实是[0,1]之间的任意数，有无穷多个的）"]},{"cell_type":"code","execution_count":8,"id":"c0ba9ea6","metadata":{"collapsed":false,"id":"1175948A07AF42CF9DC8FB5DF07D5FF4","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/plain":["<seaborn.axisgrid.FacetGrid at 0x7f2bb6dacfd0>"]},"execution_count":8,"metadata":{},"output_type":"execute_result"},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/1175948A07AF42CF9DC8FB5DF07D5FF4/sjznucr61d.png\">"],"text/plain":["<Figure size 900x900 with 9 Axes>"]},"metadata":{},"output_type":"display_data"}],"source":["# 导入必要的库\n","import numpy as np\n","import pandas as pd\n","from scipy.stats import binom\n","import seaborn as sns\n","import matplotlib.pyplot as plt\n","\n","# 设置二项分布的参数\n","n = 50  # 总试验次数\n","p_values = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]  # 不同的 p 值列表\n","k = np.arange(0, 51)                                      # 创建一个包含从0到50的整数的数组\n","\n","# 创建一个包含 'Y' 列的 DataFrame\n","dist_all_pi = pd.DataFrame({'Y': k})\n","\n","# 计算每个 'Y' 对应的概率，并将结果存储在相应列中\n","for p in p_values:\n","    column_name = f'{p}'\n","    dist_all_pi[column_name] = dist_all_pi['Y'].apply(lambda x: binom.pmf(x, n, p))\n","\n","# 使用 stack() 和 reset_index() 转换数据为长格式\n","melted_data = dist_all_pi.set_index('Y').stack().reset_index()\n","melted_data.columns = ['Y', 'p', 'prob']\n","\n","# 创建一个 FacetGrid 对象，用于绘制子图\n","plot_all_pi = sns.FacetGrid(melted_data, col='p', col_wrap=3)\n","\n","# 使用柱状图绘制概率分布\n","plot_all_pi.map(sns.barplot, 'Y', 'prob', color=\"grey\", order=None)\n","\n","# 设置 x 和 y 轴的刻度和范围\n","plot_all_pi.set(xticks=[0, 10, 20, 30, 40, 50],\n","                yticks=[0.00, 0.05, 0.10, 0.15],\n","                ylim=(0, 0.20))\n","\n","# 设置 y 轴标签\n","plot_all_pi.set_ylabels(\"$f(y|ACC)$\")\n","\n","# 设置子图的标题模板\n","plot_all_pi.set_titles(col_template=\"Bin(50,{col_name})\")\n","\n","# 显示x=30的点\n","# for ax in g.axes.flat:\n","#     ax.scatter(x=30, y=0, color='red', marker='o', s=60)\n","# g.show()"]},{"cell_type":"markdown","id":"a3622f7d","metadata":{"id":"DADE464062E64E72B663B758F432575C","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"subslide"},"tags":[]},"source":["正确率高时，正确判断次数高的情况更可能出现；正确率低时，更可能出现的是正确判断次数低的情况，这很符合我们的直觉。  \n","\n","我们可以**仅仅关注各个正确率下，正确判断次数$Y=30$发生的概率（下图黑点）**  \n","\n","![Image Name](https://cdn.kesci.com/upload/sjxndw539i.png?imageView2/0/w/960/h/960)  \n","\n","\n","\n"]},{"cell_type":"markdown","id":"6796cffe","metadata":{"id":"3ED51F7F861C48A995C5CF91A956DC89","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["在随机点运动任务中，假设参与者在 50 次试验中有 30 次正确（即  $Y = 30$ ），我们可以写出在不同正确率ACC下发生这一结果的可能性。具体公式如下：  \n","\n","$$  \n","f(Y=30|\\text{ACC}=0.1) = \\binom{50}{30} 0.1^{30} (1-0.1)^{20}  \n","$$  \n","$$  \n","f(Y=30|\\text{ACC}=0.2) = \\binom{50}{30} 0.2^{30} (1-0.2)^{20}  \n","$$  \n","$$  \n","…  \n","$$  \n","$$  \n","f(Y=30|\\text{ACC}=0.8) = \\binom{50}{30} 0.8^{30} (1-0.8)^{20}  \n","$$  \n","$$  \n","f(Y=30|\\text{ACC}=0.9) = \\binom{50}{30} 0.9^{30} (1-0.9)^{20}  \n","$$"]},{"cell_type":"markdown","id":"2f0705f5","metadata":{"id":"CC2ACD39B50A4D49A9EC92D6D4531493","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["### 似然函数  \n","\n","将这些不同正确率下  $Y = 30$  发生的相对可能性组合在一起，就构成了一个关于正确率的似然函数（$L(\\text{ACC} | Y=30)$）：  \n","\n","$$  \n","L(\\text{ACC} | Y=30) = \\binom{50}{30} \\text{ACC}^{30} (1-\\text{ACC})^{20} \\; \\; \\text{for} \\; \\text{ACC} \\in [0,1]  .  \n","$$  \n","\n","随机点运动任务中的正确率:  \n","\n","通过绘制  $Y = 30$  时，不同ACC下的相对可能性（即似然函数），我们可以观察到某个特定的ACC值最大化了这个似然函数。例如，可能在$ACC= 0.6$时，似然函数达到最大值，说明  $Y = 30$  这一结果最有可能在$ACC = 0.6$  的情况下出现。  "]},{"cell_type":"code","execution_count":9,"id":"7840b694","metadata":{"collapsed":false,"id":"FB9CB58C31F3413A8E0E66B89D174FAA","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/FB9CB58C31F3413A8E0E66B89D174FAA/sjznuctnp6.png\">"],"text/plain":["<Figure size 1000x400 with 1 Axes>"]},"metadata":{},"output_type":"display_data"}],"source":["import seaborn as sns\n","from scipy.special import comb\n","\n","# 定义似然函数\n","def likelihood(ACC, Y=30, N=50):\n","    return comb(N, Y) * (ACC**Y) * ((1-ACC)**(N-Y))\n","\n","# 定义ACC范围在[0, 1]之间\n","ACC_values = np.linspace(0, 1, 1000)\n","\n","# 计算每个ACC对应的似然值\n","likelihood_values = likelihood(ACC_values)\n","\n","# 设置Seaborn的绘图样式\n","sns.despine() # 移除图的上边框和右边框\n","\n","# 绘制似然函数，使用Seaborn的绘图功能\n","sns.lineplot(x=ACC_values, y=likelihood_values, color=\"grey\", label=\"Likelihood L(ACC | Y=30)\")\n","\n","# 设置图表标题和标签\n","plt.xlabel(\"ACC\")\n","plt.ylabel(\"Likelihood\")\n","\n","# 在ACC=0.6处画出一条虚线\n","plt.axvline(x=0.6, color='red', linestyle='--', label=\"ACC=0.6 (Max)\")\n","\n","# 显示图例\n","plt.legend()\n","# 展示图表\n","plt.show()"]},{"cell_type":"markdown","id":"908ac28a","metadata":{"id":"145A95E96DFD498B98DA637DE36D6558","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["## Part3. 后验模型"]},{"cell_type":"markdown","id":"b8d8c21f","metadata":{"id":"E9B935458F724602A8CA692EE9796BAB","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["当我们有了先验与似然两种信息，可以尝试进行推断后验：  \n","\n","$$  \n","\\begin{align*}  \n","Y | ACC & \\sim \\text{Bin}(50, ACC)  \\\\  \n","    ACC & \\sim \\text{Beta}(70, 30). \\\\  \n","\\end{align*}    \n","$$  \n","\n","* 注意，在这里，为了方便在视觉上对比先验和似然，似然函数被缩放为相加和为1"]},{"cell_type":"code","execution_count":10,"id":"6d6b3f4d","metadata":{},"outputs":[{"data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAqcAAAF0CAYAAAAep5AYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABerUlEQVR4nO3dd3hc1ZkG8Pfcqepyr8IFMDE2YMAQOoGwsIsxMVlaAiwtBAIGAwntSYBAIA4JSQgh1AWHLD0sZgkYYkPAGIwxNu69SFbv0jRNu/ec/WMkYVVrRjNz74ze3/MIpJmrez/N9cx8c+75viOUUgpERERERBagmR0AEREREVEHJqdEREREZBlMTomIiIjIMpicEhEREZFlMDklIiIiIstgckpERERElsHklIiIiIgsg8kpEREREVlGxienSil4vV5wLQEiIiKizJfxyanP50NRURF8Pp/ZoRARERHRIGV8ckpERERE2YPJKRERERFZBpNTIiIiIrIMJqdEREREZBlMTomIiIjIMuxmB0BERERdKaWg6zoMwzA7FKIBsdlssNvtEEIMel9MTomIiCwkEomgpqYGbW1tZodCFJfc3FyMGzcOTqdzUPthckpERGQRUkqUlpbCZrNh/PjxcDqdSRmJIkolpRQikQgaGhpQWlqKQw89FJqW+MxRJqdEREQWEYlEIKVESUkJcnNzzQ6HaMBycnLgcDiwb98+RCIRuN3uhPfFgigiIiKLGcyoE5FZkvXvlv/6iYiIiMgymJwSERGRKcrKyiCEwPr1680OhSyEc06JiChtoobEzoYAarwhKAWMK3ThsNH5cNiyb6xEKQVpVEPqNRBaPmyOQyBE4m+7UnoBGUxihP3QcqBphSk/TElJCWpqajBy5MiUH4syB5NTIiJKOUMqfLSrAW+sr0ZZcxBhXQIAXHYNk4bl4KKjxuOsaSNhz5IkVakIIoH3oYe+gFJBADbYHFPhyr8Imn1M3PuT0otQ6xNQ0pP8YHshtCK4i+enNEGNRCJwOp0YO3ZsUvZD2SM7XgWIiMiyWoNR/PrDnfj1h7tR2tSGsQVOTB+dh+mj8zCuwIWy5jb85qNdePjDXWgNRs0Od9CUkgj730Y0+C9AuKDZDoJmGwUjshMh398gjeb4dyqD7YmpE0IUpvQLcMaOFeco7Xe+8x3Mnz8f8+fPR3FxMUaMGIFf/OIXUEoBACZPnoyHHnoIV111FYqKinDdddf1ell/+fLlOP744+FyuTBu3Djcfffd0HW9x3Fuv/12jBw5Ev/2b/8W/+NJlsbklIiIUqYpEMEvP9iB97c3YFyBC4eMzEOeM7aKjBACuU4bDhmZh/GFbizd0YBffrADTYGI2WEPSjS0CnroC2jaCGhaUfvf6oJmL4ERLUfE/zaUSmzlJyHcEFpOar9E4i2AXnzxRdjtdnz55Zd4/PHH8cc//hH//d//3Xn/7373O8ycORNr167Fvffe2+P3q6qqcO655+K4447Dhg0b8NRTT+H555/HQw891OtxPv/8czzzzDMJx0vWxMv6RESUEoGwjoUf7cLq8lYcOioXbrutz20L3HYcOjIPq8tb8ZuPduGX5xyGPFfmvUVJownRtqUQcEJo+V3uE8IGm20s9Mgm2MLr4HDPNinK1CkpKcEf//hHCCFw2GGHYdOmTfjjH/+I6667DgBw5pln4mc/+1nn9mVlZV1+/8knn0RJSQmeeOIJCCHwrW99C9XV1bjrrrtw3333dbYqOuSQQ/Db3/42bX8XpRdHTomIKOmkVHjy8zJ8UdaCQ0b2n5h2cNk1HDIyFyvLWvDMF/s6LwdnkmjwUyijCcI2qtf7hZYDQEM0+C8oGUpvcGlwwgkndFnR6sQTT8SuXbtgGLGR4tmz+0/It23bhhNPPLHLPk4++WT4/X5UVlZ23nag/VBmY3JKRERJ98H2eizZVo+JRW7kOA6cmHbIcdgwodCNd7fWYdnOhhRGmHxSr4MeXgOhFUOIvt9eNdtoSL0aemRDGqOzhry8vH7vV0r1WK6140PK/rcfaD+U2ZicEhFRUlW0BPHC6nI4bQLFOY64f39YrgM2IfDClxWo9WbO6GI09BWU4YXQivvdTggHABuiwc8TnntqVatWrerx86GHHgqbbWAfUA4//HCsXLmyy6j5ypUrUVBQgAkTJiQ1VrIuJqdERJQ0Uio8/2U5arxhHDQsJ+H9TB6Wg/LWIP76VUVGXN5X0t8+alrQY+SvN5ptBKReBSO6Kw3RpU9FRQVuv/127NixA6+++ir+/Oc/Y8GCBQP+/RtvvBEVFRW4+eabsX37dvzf//0f7r//ftx+++1c0nUIybzZ5kREZFmflzXj071NOKg4B9oAkrS+aJrAxEI3/rWrEWdNG4XZJcXJCzIF9MgWKNkCzTZxQNvHKuJ16OH1sDu/NeDjKBUCZIJBxnOMBP3Xf/0XgsEgjj/+eNhsNtx888348Y9/PODfnzBhApYsWYI77rgDRx11FIYPH45rr70Wv/jFLxKOiTIPk1MiIkqKUNTAS2sqIRVQ6B7828uwXAdq/WG8vLYSs8YXWrZBv1IKeuhrABqEGPj8WiEKYUS2QBoeaLai/jfWciC0IijpgVKpb7UltCJAi3/k2+Fw4LHHHsNTTz3V477ulflArPdp95Hx008/HatXr+7zGJ988knccVFmYXJKRERJsWxnA7bV+TF1RG7S9jmpOAfrqrz4rLQZ3znEmktcKqMGUt8HTRsW1+8JrRDSqIQR3QnNdly/22paIdzF87Nu+VKi3jA5JSKiQQuEdby5oQYuuwaXPXkjnLnO2Ejk3zfU4OQpw+Gw4OipHtkOpdogtN7bR/UlNsoqYES2wuHuPzkFYgkqmDDSEMDklIiIBu3DXY3Y0xTAoSOT3+KnpMiNLbU+rCxtxukWGz1VSkGPbERsWdH459gKUQgjumtgl/YtjpfbKVms9xGUiIgySjBiYPHGGuTYbXCmYGQz12mDgsLbm2shpbUq96VRA6XXQNMSSyyFVgAl/ZD63iRHRpS5mJwSEdGgrChtwt7mNkwoSnxN9gOZUOjGphof1ld7UnaMRMjoHigZBERi82xjl/YV9Eh2tZQiGgwmp0RElDBDKvxjSx1sAkmda9pdgcuOUFTig23WWjVKj2wFhD2hS/odhMiDjO5ISxU+USZgckpERAlbV+XB1lofxhembtS0w5gCJ1aWNaOiJU0V6wcgDQ+kXgGhFQxqP0Ir6NwXETE5JSKiQVi6vQFhXSHflfr62hG5DjQHo/hkT2PKjzUQUi+Fkj4IkT/IPTkBRGBE9yUjLKKMZ3pyqus6fvGLX2DKlCnIycnB1KlT8eCDD0LKFC+BQUREg1LZGsQX+5oxpsCZluMJIVDosmPpjgaEouavSW9ESwEgrsb7vYlNCbDD4LxTIgAWaCX1yCOP4Omnn8aLL76IGTNmYM2aNbj66qtRVFQU13q8RESUXp/uaUJTIIoZYwc7cjhwYwtcKG1uw5qKVpwydUTajtudUgaMyHYIEf8qSr0RIh/SqIQmA31u0xyIwB/Rk3K8A8l32jE8Lz0fOoi6Mz05/eKLL/C9730Pc+bMARBbyuzVV1/FmjVrTI6MiIj6EtEllu1sQIHLBm0QxUDxctk1SAV8vLvR1ORUGnVQshVCDG6+aQeh5UMaNTD0WvR2UbM5EMGCxZvREEhP0dSoPCf+dMHMuBPUpqYmTJ8+HatXr8bkyZNTE9x+Jk+ejFtvvRW33nprwvv45S9/ibfffhvr168HAPzsZz9DJBLB448/HtexhRBYvHgx5s2bh7KyMkyZMgXr1q3DrFmzEoqrv/198sknOOOMM9DS0oLi4uKE9p+o7o9XKpienJ5yyil4+umnsXPnTkybNg0bNmzAZ599hscee8zs0IiIqA/rqjwoaw7ioOLUF0J1NyrPidXlrajzhTGmwJX24wOAjO6DUkEIbXRS9ieEHVAGpF4LYHyP+/0RHQ2BCHLsWueqWanSFjHQ0D5KG29yunDhQsydOzctiWmq3HnnnTj44INx2223YcqUKQP+vZqaGgwbFt8StgNVUlKCmpoajBxprUUoUsX05PSuu+6Cx+PBt771LdhsNhiGgYcffhg/+MEPet0+HA4jHA53/uz1etMVKhERtVu+pxFRQ8HtSG2i1JsRuQ5sa/Bj1b4WfG/m2LQfHwAMvQyANqgWUj0IR3vFfs/ktEOu05aW4rOgHn/dRzAYxPPPP48lS5akIKL0GT16NM4++2w8/fTTeOSRRwb8e2PHpu7fos1mS+n+rcb0gqjXX38dL730El555RV8/fXXePHFF/Hoo4/ixRdf7HX7hQsXoqioqPOrpKQkzRETEQ1tTYEIVpa1YGSew5Tja5qAU9Pwr12NUCr9K0YppUNG90Ak2Hi/L0LkQRr1gMrMguD3338fdrsdJ554YudtLS0tuOyyyzBq1Cjk5OTg0EMPxaJFizrvr6ysxKWXXorhw4cjLy8Ps2fPxpdffgkA2LNnD773ve9hzJgxyM/Px3HHHYcPP/yw3xg8Hg9+/OMfY/To0SgsLMSZZ56JDRs2dNnmN7/5DcaMGYOCggJce+21CIVCPfZz/vnn49VXX43r7xdC4O233+71PiklrrvuOkybNg379sW6MvzjH//AscceC7fbjalTp+KBBx6Arvc+p7isrAxCiB6X0teuXYvZs2cjNzcXJ510Enbs2NHl/qeeegoHH3wwnE4nDjvsMPzP//xPl/vLy8vxve99D/n5+SgsLMTFF1+Murq6LtsM5PFKNtOT0zvuuAN33303Lr30UhxxxBG44oorcNttt2HhwoW9bn/PPffA4/F0flVUsC8cEVE6rS5vQVMggpEmFsyMzndiR70fexrb0n5sadRCSS+EyEvqfoWWB6WCUDC/E0EiPv30U8yePbvLbffeey+2bt2K999/H9u2bcNTTz3VeWna7/fj9NNPR3V1Nd555x1s2LABd955Z2e3Hr/fj3PPPRcffvgh1q1bh3POOQdz585FeXl5r8dXSmHOnDmora3FkiVLsHbtWhxzzDH47ne/i+bmZgDAG2+8gfvvvx8PP/ww1qxZg3HjxuHJJ5/ssa/jjz8eFRUVnYnkYEQiEVx88cVYs2YNPvvsM0yaNAn//Oc/cfnll+OWW27B1q1b8cwzz+Cvf/0rHn744bj2/fOf/xy///3vsWbNGtjtdlxzzTWd9y1evBgLFizAT3/6U2zevBnXX389rr76anz88ccAYo/XvHnz0NzcjOXLl2PZsmXYs2cPLrnkks59DPTxSjbTL+u3tbVB07rmyDabrc9WUi6XCy6XOXOMiIiGOqUUPtndBJumwaalrxCqu0K3HRWeEFaXt+CQUclNEg9E6pVQKgyhJXe+rRAOAD4gQ5PTsrIyjB/fdUpCeXk5jj766M6kdf+5qK+88goaGhrw1VdfYfjw4QCAQw45pPP+o446CkcddVTnzw899BAWL16Md955B/Pnz+9x/I8//hibNm1CfX19Z57w6KOP4u2338abb76JH//4x3jsscdwzTXX4Ec/+lHnPj/88MMeo4ETJkzo/JsmTZqU6EMCv9+POXPmIBgM4pNPPkFRUREA4OGHH8bdd9+NK6+8EgAwdepU/OpXv8Kdd96J+++/f8D7f/jhh3H66acDAO6++27MmTMHoVAIbrcbjz76KK666irceOONAIDbb78dq1atwqOPPoozzjgDH374ITZu3IjS0tLOq9D/8z//gxkzZuCrr77CcccdN+DHK9lMHzmdO3cuHn74Ybz33nsoKyvD4sWL8Yc//AEXXHCB2aEREVE3la0hbK71YZRJl/Q7CCGQ67Bh+Z4mSJneS/tGtKIzhuTTAJWedlHJFgwG4XZ3Tdh/8pOf4LXXXsOsWbNw5513YuXKlZ33rV+/HkcffXRnYtpdIBDAnXfeicMPPxzFxcXIz8/H9u3b+xw5Xbt2Lfx+P0aMGIH8/PzOr9LSUuzZswcAsG3bti7TDgD0+BkAcnJiLcLa2gY3Mv+DH/wAfr8fS5cu7UxMO2J98MEHu8R53XXXoaamJq5jHnnkkZ3fjxs3DgBQX18PIPa3nnzyyV22P/nkk7Ft27bO+0tKSrpMj+x4rPffZiCPV7KZPnL65z//Gffeey9uvPFG1NfXY/z48bj++utx3333mR0aERF1s7q8Ba3BKMYXpq+3aV9G5TlR1tyGXY0BHDY6PfEoJdvnmyanv2l3QuRAwTBlLu1gjRw5Ei0tLV1u+4//+A/s27cP7733Hj788EN897vfxU033YRHH320MwHsyx133IF//vOfePTRR3HIIYcgJycHF154ISKR3ttpSSkxbtw4fPLJJz3ui7fdUsc0gFGjRsX1e92de+65eOmll7Bq1SqceeaZnbdLKfHAAw/g+9//fo/f6Z7g98fh+OZDYseHpf2vPHf/AKWU6rxt/+/72sYspo+cFhQU4LHHHsO+ffsQDAaxZ88ePPTQQ3A62fyXiMhKlFJYvqcZLruW1t6mfcl32eCPGPiqvOXAGyeJMhrb+5umZiqBgLu9ICrzRk+PPvpobN26tcfto0aNwlVXXYWXXnoJjz32GJ599lkAsVG/9evXdyaC3a1YsQJXXXUVLrjgAhxxxBEYO3YsysrK+jz+Mcccg9raWtjtdhxyyCFdvjrmuU6fPh2rVq3q8nvdfwaAzZs3w+FwYMaMGQP983v1k5/8BL/5zW9w/vnnY/ny5V1i3bFjR484DznkkB5THRM1ffp0fPbZZ11uW7lyJaZPnw4gNkpaXl7epXZn69at8Hg8ndsM9PFKNtNHTomIKDOUNrdhV6MfoyyyclDHpf3PSlvww2MmQkvDHFhpVLf3N01Rv0nhAqD6vLTfFkn9fNREj3HOOefgnnvuQUtLS2e/z/vuuw/HHnssZsyYgXA4jHfffbcz8fnBD36AX//615g3bx4WLlyIcePGYd26dRg/fjxOPPFEHHLIIXjrrbcwd+5cCCFw77339ru0+VlnnYUTTzwR8+bNwyOPPILDDjsM1dXVWLJkCebNm4fZs2djwYIFuPLKKzF79myccsopePnll7FlyxZMnTq1y75WrFiBU0899YCjuwNx8803wzAMnHfeeXj//fdxyimn4L777sN5552HkpISXHTRRdA0DRs3bsSmTZvw0EMPDfqYQGzk+eKLL+4sCvvHP/6Bt956q7PjwVlnnYUjjzwSl112GR577DHouo4bb7wRp59+eucc4YE+Xslm+sgpERFlhrUVHnhDOorc1hnXGNl+aX9vU3qq9g29GlAKQqTm7fObS65dk9N8px2j8pwI6hJNbdGUfgV1iVF5TuQ74zvPRxxxBGbPno033nij8zan04l77rkHRx55JE477TTYbDa89tprnfctXboUo0ePxrnnnosjjjgCv/nNb2CzxXrn/vGPf8SwYcNw0kknYe7cuTjnnHNwzDHH9PvYLVmyBKeddhquueYaTJs2DZdeeinKysowZswYAMAll1yC++67D3fddReOPfZY7Nu3Dz/5yU967OvVV1/FddddF9ff359bb70VDzzwAM4991ysXLkS55xzDt59910sW7YMxx13HE444QT84Q9/GFTxVXfz5s3Dn/70J/zud7/DjBkz8Mwzz2DRokX4zne+A+Cb1lfDhg3DaaedhrPOOgtTp07F66+/3rmPgT5eySZUJk5s2Y/X60VRURE8Hg8KCwvNDoeIKCsppbBg8RZsqfXikJHprY7vj1IKW+r8uPW0qbhoVt/N65Ml2PoXGNFS2OypOVY4koOqhqMwZfIk5OZP6HJfc/uqTemQ77THvToUACxZsgQ/+9nPsHnz5qRdnk639957D3fccQc2btwIu906H8QyQSgUQmlpKaZMmRLX3Nnu+KgTEdEBlbcEsbvRb2pv094IIeCya1hZ1owLjxqX0kIOJdsgjbqkN9/vSUBBh1Kyywjt8DxnQgljOp177rnYtWsXqqqqMnaRnEAggEWLFjExNREfeSIiOqB1VbFL+hOLktvbMxlG5DqwsyGAak8IE4pTU0UPANKogZIBaLbBVXAfkBCxoigVbZ+DmlkWLFhgdgiDcvHFF5sdwpCXmWPuRESUVitLW+CwCdNbzPSmOMcBb0jHxhpvSo8j9VrEquhT3eNVAFBQGVixT5QMTE6JiKhftd4QttX7MCLXmpeUNSGgCWBNRWtKj2Po1QBS1Xz/G517V9GUHofIqpicEhFRvzZUe9EajGJYjrmrQvWnOMeBdZVe+EKpGW1USkHqpRAiHdMaFJQSgOq92TyRVSWrxp7JKRER9Wt1eQuEEGnpI5qo4bkONLVFsLk2NZf2lfRASQ+Q4mIouz0CpSSCIQMKUSjVd19PIqvpWHp1/5WrEsGCKCIi6pMvpGNdpdfSo6YA4LRp0KXCxmovTpzc+1rtgyGNWijZBs1WdOCNB8GmGSjIrUZDgx3AMOTl+6Fp1pxOQdRBKYW2tjbU19ejuLi4s1dtopicEhFRnzbVeNHUFsHUEalunzR4BS47Vu1rwXUnTEr6KK8y6gFICJH6t82RxeUAFOrrgxCNbRCCySllhuLiYowdO3bQ+2FySkREfdpQ7YWhFJw2688CG57jQLU3jD1NARw6Kj+p+zb0KkClZ1qDEMCoYRUozlsFOP8dzpxT03JcosFwOByDHjHtwOSUiIh6pRsSX5Q1oyDOZSzNku+yodwTxJZaX1KTU6UUZHQfhJbeHq+apmDXKga10g5RJrL+R2EiIjLFrsYAan1hDM+19nzTDkII2ITA15WepO5XyVYo6YUQqWvw3xsh3JBGNZQy0npcIrMxOSUiol5trvEhEDGQ50zOpbp0GJbjwKYaX1JbSkmjDkoFgbQnpzlQMgAlW9J6XCKzMTklIqJerS5vhdOiq0L1ZViOAy3BCLbV+ZK2z3QWQ3Uh3FAqCGk0pPe4RCZjckpERD00BSLYUe/HcIuuCtUXp12DbihsqU1ecmrotdhv3aa0EcIGQEEZjWk/NpGZmJwSEVEPm2u8aAlGUWzx/qa9yXHY8FVFa1JWq4mtDFUOIVxJiCwxUq837dhEZmBySkREPWys8UIpwG7hVaH6UpxjR3lLELW+8OB3pvxQsjXt8007CLgg9QpTjk1kFianRETUhW5IfFXeigJX5hRC7a/I7YAnpGNbnX/Q+5JGPZQMpr1Sv5NwQ8pmKBky5/hEJmBySkREXextasuoFlLd2TQBBYXNNd5B7ytWjKTDrLbgQrgBFYKUnHdKQweTUyIi6mJrnQ9t0cxqIdVdgdOOtZUeSDm4eadSrwOgzOtYIJxQKsyiKBpSmJwSEVEX66q8sInMaiHVXXGOA3W+MMpa2ga1H6lXAjCvY4EQsbdpaTSZFgNRujE5JSKiTr6Qjk3VXgzLwCr9/RW4bPCHdeyoT3zeqVIRSKPevPmmnQSkUWdyDETpw+SUiIg6ba/3Z2wLqf0JIaAJgc01ifc7lUYjoIKxeZ8mEsIFqVclpTUWUSZgckpERJ221vkQNSRc9sx/eyhw27GuygPdkAn9vjIaoFQYMLHHKYDYSlHSA6jBTVEgyhSZ/+pDRERJoZTC2opWuB3Z8dZQ7LajMRDB3qbEkrqOZUM75n2aJVaxH+a8UxoysuMViIiIBq3BH0Fpc1vGzzftkOe0oS1qYEdDYvNOpV4DM5Yt7ckJpSJQkskpDQ1MTomICACwrc6H1qCOQnd2JKcd80631MY/7zS2bGmVBYqh0Nk1gSOnNFQwOSUiIgDA1jo/lFIZuWRpXwpddqyv8iIa57xTJb1QygeYXAz1DQFp1JsdBFFaMDklIiIopbCmojWjG+/3pshtR3NbBKVxzjtVshFKhkyv1O8Qq9ivNjsMorRgckpERKhsDaHGG8r4FlLdxeadSuxsCMT1e9JohJnLlvbQXrGvJCv2KfsxOSUiImyr98Eb1lHotkgyliRCCAjE5tPGQ+odlfrWmOIg4AJUCFI2mx0KUcoxOSUiImxtLxrSLJKMJVOh24711fH1O5V6FSwzagoAor1i32ByStmPySkR0RBnSIWvKz0odFkoGUuiIrcdjYEo9rUEB7S9UgakrLPMfFPgm16rkskpDQFMTomIhriy5jbU+yMoypIWUt3lO20IRAzsGuC8UyWbAdlmqeQUAKDAin0aEpicEhENcTvq/fCHdRS4sqtSv0PnvNP6gc07lUYTlApZqI1UjBAuKKPG7DCIUo7JKRHRELel1gchrFP8kwr5ThvWV3khpTrgtspoBKAghMWSdeGCNJqhVMTsSIhSiskpEdEQphsS66uzd75phyK3HfX+MKo8oQNuK42GNEQUPyFiFfvKaDE7FKKUYnJKRDSElTa3odEfybr+pt0VuOzwh3XsajzwvFOpV0LAmYao4iRcUCrMdlKU9ZicEhENYTsbAghEZdatDNWdpgkoBexs8Pe7nVKRWAN+i803BRCbZqAUR04p6zE5JSIawrbW+iCQ3fNNO+S2zztVqu95p8poBlQodgndioSCNJrMjoIopZicEhENUbohsa7Kk3WrQvWlyG1HlSeIxkDfBUWxSv2wJUdOYxyQRp3ZQRClFJNTIqIhqrS5DU2BCIqHSHJa6LLDE9L77XeqZEelvjXfHmPtpGr7Hf0lynTWfPYREVHKxeabGlk/37SD3aZBSoXd/RRFGXoDYOW8T7igVBBKes2OhChlmJwSEQ1RsfmmYkjMN+3gtGvYVNN3M35lVFl3vinaR05lGEqyKIqyF5NTIqIhKNbf1IvCLF0Vqi9FLjt2NwYQCOs97lMqHCs2snByCjgARJmcUlZjckpENATta4kVBhW5s7u/aXeFbju8oSj2NLX1uE8ZTYAKQ1i2GOqbrgqS7aQoizE5JSIagnY2+BGIGMgfYiOnLruGsCF7nXcaWxo0ZPGRUwDKuqtYESUDk1MioiFoe71/yPQ33Z8QApoQ2FHfsxm/kk3t21j7rbGjYp8oW1n7GUhEREknpcL6Ki/yh0iVfncFLjs21XihG7LL7Zav1O8gXO2jvH33ayXKZExOiYiGmIrWIOr9YRTnDI3+pt0VuexobouiojXU5XarV+p3EMIFqDCU9JgdClFKMDklIhpidjUG4A8byHcNzeQ0z2VDIGJgT9M3806VimRApX474YRSYSgWRVGWYnJKRDTEbK/zQykFbYjNN+2gCQGlgJ31+yWnGVCp30EIOwAJyXZSlKUskZxWVVXh8ssvx4gRI5Cbm4tZs2Zh7dq1ZodFRJR1lFJYX+0ZclX63eU6NWyq8XYuAxqbwxnOjJFToL1iv9nsKIhSwvRrOi0tLTj55JNxxhln4P3338fo0aOxZ88eFBcXmx0aEVHWqfGGUeMNo3CI9TftrtBlR6UniOa2KEbkOdsr9ZXlK/U7CQ2K7aQoS5menD7yyCMoKSnBokWLOm+bPHmyeQEREWWxXQ1+eEM6xhVkyAhhihS47WhsasPepjaMyHNC6pmV6AnhgjRqoZQacu3AKPuZ/hHxnXfewezZs3HRRRdh9OjROProo/Hcc8/1uX04HIbX6+3yRUREA7OzIQClFGza0E5onDYNulKdRVHSqIaA0+So4iBcUNIHqNCBtyXKMKYnp3v37sVTTz2FQw89FP/85z9xww034JZbbsHf/va3XrdfuHAhioqKOr9KSkrSHDERUWZSSmFdlQc5jqE937SDXQhsqfW1V+o3Zs58UyCWSKswi6IoKwnVMRvcJE6nE7Nnz8bKlSs7b7vlllvw1Vdf4YsvvuixfTgcRjgc7vzZ6/WipKQEHo8HhYWFaYmZiCgTNfrDuPq19XDZNIzIy6BRwhSp8Ybgdtjw4iUTIP2PQYgCCC3X7LAGRCkD0qiEu/DHsLtmmB0OUVKZPnI6btw4HH744V1umz59OsrLy3vd3uVyobCwsMsXEREd2O7GNnhDOgrdppcbWEKBy47WYBSlTfWZVakPQIjY6LeSreYGQpQCpienJ598Mnbs2NHltp07d2LSpEkmRURElJ12NvhhSAWHzfSXfkvIc9oQjErsbWxGrFI/86Y7sJ0UZSPTX6Fuu+02rFq1Cr/+9a+xe/duvPLKK3j22Wdx0003mR0aEVFW2Vjthctu+su+ZQghoJTC7oZMLay1QRn1ZgdBlHSmv0odd9xxWLx4MV599VXMnDkTv/rVr/DYY4/hsssuMzs0IqKs4Qvp2N0U4CX9bvKcNmysaYNSmdf3NdZOqh4ml44QJZ0lXqXOO+88nHfeeWaHQUSUtXY3BuAN6TioOMfsUCylwG1Djc9AUzAPYzItPxUuKBUAVAAQ+WZHQ5Q0po+cEhFR6u1qDCCiS17W76bAKeELA/s8mVdcK0R7OymD7aQou/BViohoCNhS6xvyjfd7Y9ciMJRCmSfP7FAS4IRSEVbsU9ZhckpElOXCuoGttT4Ucb5pTyoMu5DY0eg2O5K4CRF7C2dyStmGySkRUZbb29SG1mAURe5Mm1SZekqFUODSsb3RiYhhdjSJUGwnRVmHySkRUZbb3RhAMGogx8GX/B5kG/KdEt6whvLWTEzeHZBGndlBECUVX6mIiLLc9jo/hBAQgnNO96cAKBVEnkOgLSpQloHJqRAuKKOB7aQoqzA5JSLKYoZU2FDtRaEr81Y/SjkVgVI6NE0DILC3JfOSUwgnlGoDlN/sSIiShskpEVEWq2gNojEQQSHnm/akQgAMAHbkOCQ21zuRaQOQQrgAFYE0Ws0OhShpmJwSEWWxXQ0B+CM68jly2oNSYQASgECBU6Laa0dLKNPeFh3t7aTY65SyR6Y9C4mIKA47G/yAEtA437SHWHIKQAgUuCR8ES3j5p2ynRRlIyanRERZSimFdVVe5Dn5Ut8r1QYglrQ7bYAugbJMnHcKsJ0UZRW+YhERZal6fwS13hD7m/ZCAVCyDcA30x00AexsysTHygGps50UZQ8mp0REWWpXgx/esI4CrgzViygUooD4JjktcEpsq3dBlyaGlQAhnFCS7aQoezA5JSLKUrsaApBKwa5xvmkPMgwoA2K/kdMCl0RLSEOVN8OSeeGCUkEo5TM7EqKkYHJKRJSlNlR74bLxZb43Ch2V+t88PnlOlZHN+GMjpxEotpOiLMFXLSKiLNQajKK0uY3zTfugZCj2zX5dDDQRm4ta2pxpj5kTQIQV+5Q1mJwSEWWhPY0BeEM6ijjftHcqhI5K/f257MDmBmf64xmEjmVp2euUskXCyWkkEklmHERElES7GwOISgmnnWMQvVEqgN7eAgucEhWtDvjCmTdPVxpMTik7JPyqNWHCBNxzzz0oLy9PZjxERJQEm2t8sLPxfq8UDCgVwf5tpDoUuCS8GdiMn+2kKJsknJzOnTsXjz/+OA4++GBccMEF+Oijj5IZFxERJSgUNbCt3odCXtLvnQrFKvVFz+TUZVOIGMi45JTtpCibJJycvvDCC6isrMTDDz+MDRs24Oyzz8b06dPxxBNPwOdjOwsiIrPsbWpDa1BnMVQfYsuWGuht5FSI2EzU3ZlWFMV2UpRFBjUZadiwYbjzzjuxZ88eLF68GCUlJViwYAEmTJiA+fPnY/v27cmKk4iIBmh3YwAh3UCOg/NNe6NkOPZNH9Me8pwKW+pckBk0CMl2UpRNkvLKJYTA+eefj0ceeQSnn346/H4/nnzyScyYMQP/+Z//ifr6+mQchoiIBmB7vQ+aEJ1V3NSNCvZ7d6FLoiloQ62/58iqdbGdFGWPQSenuq7j1VdfxSmnnILZs2dj7969eOSRR1BWVobHHnsMK1aswH/9138lI1YiIjoAQypsqPaiwJlJiVV6KdWG3i7pd8hzSvjCGvZl0LzTb9pJtZobCFESJDxbvqqqCs888wyee+451NXV4dRTT8Ubb7yBCy64AJoWy3lvvvlmTJgwAZdffnnSAiYior6VtwTRFIhiWE7mJFbppCDb55z2nZzaNUAphdIWB04sCaUvuCSQRrPZIRANWsLJ6eTJk2G323HppZdiwYIFmDVrVq/bTZ06FWPGjEn0MEREFIddjX4EIjpKit1mh2JNKgwoHUL032jfYQO2Z1gzfsABaXAaHWW+hJPT+++/H9dffz1GjRrV73azZs1CaWlpoochIqI47KwPQAHQON+0V7FRU4kDzWorcEnsanYgpAu47ZlRGSWEE8qItZPifGPKZAnPOT3ooIM6L99319zcjL/97W8JB0VERPFTSmF9tQf5nG/aNxUCoABxgOS0fd5peWsG9YoVLijVxnZSlPESTk6vvvpq7Nmzp9f7SktLcfXVVyccFBERxa/WF0atN4xC9jftk1IDm0Oa41AIRgVKM6ooiu2kKDsknJz2twpFKBSCzcZP7kRE6bS7MQBvSEehK4NG+9JMyf4r9TsIEfvP3oxqxu8A20lRNojrFay8vBxlZWWdP69btw6hUNdPocFgEM8++ywOOuigpARIREQDs7PeD6kUbBrnG/ZGQUGpIAY6LpPrkNhc74JSffbrtxTRPlWBySlluriS00WLFuGBBx6AaG/ufOONN/bYpmNE9U9/+lNyIiQiogHZUO2Fm6tC9U2FAWVAiIG99RW4JGr9NjQFNYzMlSkOLnnYTooyXVzJ6cUXX4yZM2dCKYWLL74Yv/71r3HooYd22cblcmHmzJmYPHlyMuMkIqJ+tLRFUNYSRBHnm/YpNt/UAOAa0PYFTonGgANlrQ6MzA2nNLbksUPqdWYHQTQocSWn06dPx/Tp0wHERlHPO+88jBgxIiWBERHRwO1ubIM3FMWU4blmh2JdKoyBVOp3cNgAQwH7WhyYPT4zklMhXFCyke2kKKMlPGv+yiuvTGYcREQ0CLsa/NANBaeNl/X7EutxGh+bALY3ZlAzfuGMLc+q/IAoMDsaooTElZw++OCD+NGPfoTx48fjwQcf7HdbIQTuvffeQQVHREQDs7HGCwcT034pGUC8TWoKXBI7Gp3QZWxZU6uLjZwGII0W2DQmp5SZhOqvJ1Q3mqZh1apVOP744/tswN+5YyFgGMagAzwQr9eLoqIieDweFBYWpvx4RERWEwjruPzldVBKYUzBwOZTDjUKCnp4XawgShv40q7+sEBT0IbHz63HlGF6CiNMDqUUpFEOd+G1sLuONDscooTENXIqpez1eyIiMs+epjZ4QlGUFA086RpyVARQ+oAr9TvkOhUqvBpKWxwZkZx2zDNlOynKZBlwkYKIiPqzq8GPiC7hyoTrzmbprNSPb4GYWMtYhdKWzOqCwHZSlMkSfiULhULwer1dbnvjjTdw991348MPPxx0YERENDBban2waYLV2f1QcVbq789tV9hUn0nTJRxsJ0UZLeHk9IorrsAtt9zS+fPjjz+OSy+9FL/97W9xzjnnYMmSJUkJkIiI+hbRJTbX+lDo5pKl/Yn1OE1MgUuiymtHaygzRqaFcELJBijF6XeUmRJ+pq1evRr//u//3vnz448/jssvvxytra34/ve/j0cffTQpARIRUd/2NgXQEoyiiMlpv5RsQ6JveQVOCV9YQ1mmXNoXLigVhJI+syMhSkjCyWlDQwMmTJgAACgtLcXevXtx8803o7CwENdeey02b96ctCCJiKh3uxoDCEYM5Drim0s5lCgg1vszzvmmHVx2ICoFylozIzmNjZxGWBRFGSvh5DQ3NxcejwcAsGLFCuTn52P27NkAALfbDb/fn5wIiYioT9vqfBCC8037pcLtlfqJJ/CaUNjVlBnJKeAEwOSUMlfC14GOOOII/OUvf8GkSZPw5JNP4owzzuh8cSwvL8fYsWOTFiQREfWkGxIbqr0odHHUtF+dlfqJFzXlOyW21DthSMDqax2wnRRluoST03vvvRfnnXceZs2aBafT2aVC/7333sMxxxyTlACJiKh3+1qCaAxEMSInU0b0zNFZqY/ER5cLXBItIRuqvHYcVGz9fqcA20lR5ko4OT3zzDOxbds2rF27FrNmzcLUqVO73Ddr1qxkxEdERH3Y1RBAIGJgUjGb7/ens1J/EFMf8p0K1V6B0lZHhiSnDki91uwgiBIyqPLOSZMmYdKkST1uv/766wezWyIiGoDt9bFqbM437Z+SAQx2zRlNxMZeS5sdOH1yMClxpZIQLkjZCKUkRAK9XYnMNOjeI/X19di3bx+CwZ5P1tNOO22wuyciol5IqbCuyosCJ+eb9kdBQakgEq3U35/LDmxucA4+qDSIVewHoaQXwlZsdjhEcUk4Oa2pqcEVV1yBjz/+uMd9SikIIWAYxqCCIyKi3lW0BlHvD6Moh/1N+6Ui7ZX6g3+cCp0S5a0OeMMaCl0Wb3AvXFDSFyuKYnJKGSbhZ+v8+fOxbt06PPLIIzjyyCPhcmXS0m5ERJltV2MA/rCBiUWcb9qvJFTqdyhwGaj0OFDaYsdRYyOD3l9qOQBEoWQLgMkmx0IUn4ST0+XLl+PRRx/F1Vdfncx4iIhoALbX+aGgoHG+ab86K/WTMO/SZQcihkBpi8PyyakQsUmy0mgxOxSiuCX8bBVCoKSkJJmxEBHRACil8HWVB/lOXtI/kM5K/STRNIWdTZkx7xRgOynKTAknpxdddBHefffdZMaChQsXQgiBW2+9Nan7JSLKJlWeEGq9YRS5mZweSDIq9feX75TYWu+EbvEpp0B7UZTBdlKUeRJ+Zbv44otx3XXXQUqJuXPnYsSIET22iacR/1dffYVnn30WRx55ZKIhERENCTsbAvCFo5hQyLn+/UlmpX6HApdESzDWjH+S1fudChek0QSljEEt3UqUboNqwg8ATzzxBP7yl790uS/ean2/34/LLrsMzz33HB566KFEQyIiGhK21/mhlICmcb5pv1Q4aZX6HfKdClVeDXtbHJZPToVwQqkglGyFsPUcQCKyqoSfsYsWLUpaEDfddBPmzJmDs84664DJaTgcRjgc7vzZ6/UmLQ4iIquLzTdtRZ6TjdUPRCWxUr9D7POAQmmzA2dMsXgzfuGCkp72dlJMTilzJJycXnnllUkJ4LXXXsPXX3+Nr776akDbL1y4EA888EBSjk1ElGmqPSHUeMMoznGYHYr1qRCSVam/vxy7wsa6TJhSYQegQxotsPGfC2WQpDxjd+zYgc8//xyBQCCu36uoqMCCBQvw0ksvwe0eWK++e+65Bx6Pp/OroqIikZCJiDLSzoYAvCEdBS4WQx1Isiv1OxS4JKp8drQErT163bGsbazXKVHmGNQz629/+xsmTpyIww8/HKeddhp27NgBIFYs9dxzzx3w99euXYv6+noce+yxsNvtsNvtWL58OR5//HHY7fZe56y6XC4UFhZ2+SIiGiq2tfc3tXG+6QHFKvWTXwhU4JLwhjXsacmE4UgBaTSYHQRRXBJOTv/+97/jqquuwjHHHIMnnngCSqnO+4455hi88cYbB9zHd7/7XWzatAnr16/v/Jo9ezYuu+wyrF+/HjYbqwuJiDp0zjd18LXxQBRk+8hp8kc3nTZAl0BZBiSnAi5Ine2kKLMkfF1o4cKFuPrqq/H888/DMAzcdNNNnfdNnz4df/7znw+4j4KCAsycObPLbXl5eRgxYkSP24mIhroqTwjVHs43HRAVaq/UT03DfJsGbGvIgGb8wgklW6FUJGWPBVGyJfyRctu2bbj00kt7vW/48OFoampKOCgiIuppR70fvnAUhZxvekDfVOqnZpS50CmxvdGJsG7t6RVCuAAVhuIyppRBEn6Fy83Nhcfj6fW+qqoqDBs2LKH9fvLJJ4mGRESU1bbXs7/pQCnZXgwlUvNYFbok6gM2lLXacdjIaEqOkRTCBSUbIWULNIwxOxqiAUl45PTkk0/uMde0w1//+ld85zvfGUxcRES0H6UU1lZ6kO/ifNMBUUEAqUvicxwKwWisGb+VxVaGUlBGs9mhEA1YwiOn9913H0455RQcf/zx+OEPfwghBN566y3cf//9+PTTT7F69epkxklENKRVtoZQ6w2h2M1L+geiACgVQCqKoToIAQihsLvJARyassMkjeRlfcogCT9zZ8+ejffffx9+vx8//elPoZTCr3/9a+zcuRNLlixhQRMRURJtr/fDFzZQwOR0AHQoFQaSuGxpb/KcChvr3JA9LyBajA3KqDM7CKIBG9Qz94wzzsC2bduwZ88e1NXVYeTIkZg2bVqyYiMionZbar0AFLQUzaHMKjIEKANCpPaSe6FLoiFgQ7XPjomFekqPNRhCuCCNOiilOhvzE1lZQslpQ0MDnnnmGXz66aeorq4GAIwfPx5nnHEGfvzjH2PECK7hS0SULIZU+LrSg3wnR00HIlapL5HKy/oAkO+UqPHZUNrisHRyCuGKTXNQfkAUmB0N0QHF/Ur30Ucf4T//8z/h9Xphs9kwcuRIKKWwY8cOfPjhh3j00UexePFinHbaaamIl4hoyNnX0oZ6fwTD2N90QJQKxr5J8SihTYvNb93d5MCpk4IpPdZgxEZOA5BGC2wak1Oyvrg+VjY0NOCSSy5BUVER3njjDXg8HtTU1KC2thYejwevvfYa8vLycOGFF7LPKRFRkmyv88Mf1lHASv2BkW1IZaX+/tx2hY11Vm9u7wQQgZKs2KfMEFdy2rEa1Oeff44LL7wQubm5nffl5ubi4osvxmeffYZoNIrnn38+6cESEQ1Fm2u9EEJwvuAAKChIFcAgSyoGrNAlUeFxoDmY2ikEgyGEABQr9ilzxPVsWrp0Ka655hpMnDixz20OOuggXH311fjggw8GHRwR0VAXNSS+rvSiiFX6A6PC7cuWpmeUudAl4Q1r2N1s/SkX0mg0OwSiAYkrOd22bRtOOeWUA2536qmnYtu2bQkHRUREMbsbA2gKRFCcw+R0IFK9bGl3DhtgKIG9zda+tC+EE8qoMTsMogGJKzltbW3F6NGjD7jd6NGj0drammhMRETUbludH8GogVwH55sOiAoBUIBI32V2h6awpd7aySmEC9JoglIWXmqVqF1cz95wOAyH48CXLux2OyKRSMJBERFRzIZqL2wa55sOlJJtaT9moVtiZ5MDgYh1z5EQbkCFoTjvlDJA3NeJduzYAbu9/1/bvn17wgEREVFMIKxjU40XxWwhNWCxZUvTO8pc6JKo8tqxp9mBI8dadGBGuKBkE6RshoYDXwElMlPcyelVV111wG24CgUR0eDtaPCjpS2KScNyzA4lI6g0LVvancumEDEE9rY4LZucxgrEJJTBdlJkfXE9gxctWpSqOIiIqJsttT5EDAmX3bptiixFhtor9d1pPawQgE3E5p3Om57WQ8dNMjmlDBBXcnrllVemKg4iItqPUgprKjxwMzEdsNjKUKlftrQ3hS6JrQ1OhHUBl12l/fgDY4dkxT5lAL7qERFZUIM/gr1NAS5ZGod0LVvam0K3RGtIQ2mLdVt+CeGCMuqglDQ7FKJ+MTklIrKgbXU+tAZ1FDE5HTAlAzDrbS3HrhCKCuyycL9TIVxQsg1Kes0OhahfTE6JiCxoc60PAGDXWFw6EAoSSrUh3ZX6HYSILRO6vcG6ySmEG0qFoWST2ZEQ9YvJKRGRxRhSYXV5K/KdbLw/YKqjGMq8y+qFbgOb6lyIGqaFcAB2ADqLosjymJwSEVlMaVMban1hzjeNQ2y+afqWLe1NkUuiOaihrNWa562jxSPbSZHVMTklIrKYLXU+BMI68l0cOR0oJUOxb0zssZ3rUGiLatjZZOFL+9AgjVqzgyDqF5NTIiKLWVfl4ZKl8ZJ+AOY+XkIAmqawzcLzToVwQ+o1UMqq7a6ImJwSEVmKL6RjYzWXLI2HgoJUASSw6GHSFTolNtZaeN6pcEEpP6D8ZkdC1Ccmp0REFrKtzofmtgiGMzkduM5iKPOnQRS5rT7v1AUlQ5AGK/bJupicEhFZyOZaH6KGgpMrQw1YrBhKh5nFUB1yHQpBS887dQKIsiiKLI2vfkREFqGUwqp9LchjC6m4KNmxMpT5b2mxfqfA1nprJqcd85gle52ShZn/TCYiIgBAeUsQFa1BtpCKlwrA7GKo/RW5DWyscyFi1Xmn0CB1VuyTdTE5JSKyiE01PvjCOgrd5hf2ZAoF1b5sqXUes45+p3ssupSpEC5Io5oV+2RZTE6JiCzi68pWaEJAYwupgVNhKBW1RDFUhxyHQjAqsLPJoiPgwg0lfazYJ8tickpEZAG+kI711V5e0o+TUm2wSjFUByEAmwZsrnOZHUqvhHCzYp8sjckpEZEFbKn1oTnAFlLxslIx1P6K3BKb651oi1pxFNwBqCik0Wh2IES9stazmYhoiNpQ7YGu2EIqbhZYGao3xW6J1pANuy14aV8IAQhAceSULIqvgkREJjOkwsqyFuQ7rVPUkwmstDJUd267QkQHdjRasygKsEEaNWYHQdQrJqdERCbb3RhAjTeE4bnWG2WzNBUEVBRCWC85BQCHDVhXY915p1KvYsU+WRKTUyIik22s9iIQMZDP5vtxia0MZcBKxVD7K3ZL7Gp2oiVovbfa2DKmASjZanYoRD1Y7xlDRDSEKKXwRVkzXDatc/UeGhgl22LfWPRxK3IbaA1p1ry0L9xQKgTFoiiyICanREQmqvWFsbMhgBF5vKQfLyV9sPLbmMMGSCWwvcGCySnsAHRW7JMlWfdZTUQ0BGyo8qI1FEWxm8lpPBSM9h6n1pxv2iHHIbGm2g2rTe2MjdILJqdkSUxOiYhM9FVFS2xVKM2al6YtS7YBSrdsMVSHYW6JSq8dVV4rxumA1CvNDoKoByanREQm8YaiWFfpZeP9BMRGTSWs/jZW4JLwhTVss+C8UyHckEYdlIqYHQpRF9Z+VhMRZbEN1V40tkXYQioBSgZi31i0GKqDJgAhFDZZcClTIdyACkIZzWaHQtQFk1MiIpOsqWiFVAoOG1+K46EAKOWD1eebdihyS6yrcSGsWyyRbq/Yl0aD2ZEQdcFXRCIiEwQjBr4oa8EwFkLFT4WhZNjy8007FLslmto07Gi01rkWQgMUmJyS5TA5JSIywcYaL+r9EYzIs95cRKuLzTe1bvP97tx2hZCuYZsVW0oJAWnUmh0FURdMTomITLC2woOoIeGy82U4XrH5pgoQmfHYCQG47ApfVVmxpZQbMlrOZUzJUjLjmU1ElEXCuoHPSptQ5M6My9JWY/Xm+70ZnmNgb4sDdQGLjfaKHCjlg5IesyMh6pRZz24ioiywucaHWl8YI3lJP26x5vsBZEoxVIdCl4Q3bMPWemudcyHcUDIIZdSbHQpRJyanRERptqaiFWFdIsdhsVG0TCADGdF8v7tYQwaFDbXWaiklhANQBouiyFKYnBIRpVFYN/Dp3mYUujIrubKKTGm+35sit8TaajdClmspBUidRVFkHZn37CYiymCbanyo8YYwKt9al3czRWy+qbB88/3eDMuRaGqzWa5qX8AFQy83OwyiTkxOiYjSaPW+VkR4ST8hChJSZk7z/e7cdoWIAWyqs1ZyCi0HSjZDSb/ZkRABYHJKRJQ2oSir9AdFBgEVzbj5pvvLcSqsqsiBtFDnJiFyoGQQkkVRZBGmJ6cLFy7Ecccdh4KCAowePRrz5s3Djh07zA6LiCjpNlR7Ue0NYyQv6SckVqWfOc33ezM8x0Cl146yFiutFuUAEGVySpZhenK6fPly3HTTTVi1ahWWLVsGXddx9tlnIxAImB0aEVFSfVHWDN2QcNszN7kyUybPN+1Q4FTwRwQ2WujSvmh/PFkURVZh+rWRDz74oMvPixYtwujRo7F27VqcdtppJkVFRJRcvpCOz0qbMSzHSiNmmUNBQUkvLPC2NShCAA4b8GWlG/OmW2cQRsAFqZeZHQYRAAs+yz2e2CoVw4cP7/X+cDiMcDjc+bPX601LXEREg7G2shX1/ggOHpFrdiiZSbVBqQiEsM6IY6KG5xjY2ehErd+GsfmG2eHEaDmQRhOUDEBoeWZHQ0Oc6Zf196eUwu23345TTjkFM2fO7HWbhQsXoqioqPOrpKQkzVESEcVvxd4mKAU4bZZ62c0YSmb+fNMOw3IkWsM2bLRQQ/5YUVQbpMFL+2Q+S71Kzp8/Hxs3bsSrr77a5zb33HMPPB5P51dFRUUaIyQiil+9L4yvKjwYmcdL+onKhvmmHTQBaELhqyq32aHsxwFAh9TrzA6EyDqX9W+++Wa88847+PTTTzFx4sQ+t3O5XHC5rPNpk4joQL4sb0FzIILpo/PNDiUjZct80/0Ny5HYUOtCS1DDsBxpdjjtRVEChl4NfoQis5k+cqqUwvz58/HWW2/hX//6F6ZMmWJ2SERESaOUwkc7G+GwCWha5o/6mUIG2uebZk9yOjzHQHNQw3pLXdp3QeqlUMpCTVhpSDI9Ob3pppvw0ksv4ZVXXkFBQQFqa2tRW1uLYDBodmhERIO2qyGA7fV+jCmwThKSaZTyI1vmm3awt7/7rqm20KV9kQtltEDJVrMjoSHO9OT0qaeegsfjwXe+8x2MGzeu8+v11183OzQiokH7rLQZvrCOQlf2jPqlmzI8yJb5pvsbliOxpsoFT8j0t2IA7UVRikVRZD7TXy15+YCIslUoauCjXY0odNk7G51TfBQMSOWHBd6ukm54joE9zQ6sr3Xh9MnmXy0Uwg4oBanXAM7pZodDQ5g1Pq4REWWhrypaUdka5CX9QVDSD2RJf9PuHO2zFFZXWunSvgapl5sdBQ1xTE6JiFLko52NkApw2flSm6hYCykFiOx8DIfnSKypdqMlaI2/T4hcGHo5lIqaHQoNYdZ4NhARZZmKliDWVLRidH72jfiliwKgZAuyqRCqu+G5BpraNKyrscbouhC5gPRD6px3SuZhckpElALL9zSiJRjFiFx2jUyYCkHJYFa1kOrOrgEQwMqKHLNDiRFuKBWGNKrNjoSGMCanRERJFooaWLqjAflOGwuhBiF2SV8Hsrwt/MhcA+trXKjzmz9CHPv3qiD1KrNDoSGMySkRUZKtLm/FvpYgxhVaqNAlAymZnS2kuhueI9EStFmm56kQOTCiu9lNh0zD5JSIKImUUnh/Wz0UWAg1GAoGpPQgG1tIdacJwG5TWF6WAyvkg0LktTfjbzI7FBqi+MpJRJREuxoCWFflwVi2jxoUJX1Z20KqN6PzDOxocKK0xQJTGERurBk/L+2TSZicEhEl0bKdDfCFdRS7s3/EL5WU9CKbW0h1V+iS8EU0rLJAz1MhNEApGNEKs0OhIWpoPOuJiNKgKRDBRzsbMTzXwUKoQVBQUEYzsrmFVHdCAPlOiX/tzUXEMDsaQAgXjOguzjslUzA5JSJKko92NaDOH8aYfF7SHxQZgFLhIXNJv8OoPAMVXjvW11rg34+WD2U0tPeZJUovJqdEREkQjBh4b2s98pw22DSOmg5GrBDKwFAaOQWAHIeCLgU+Lcs1OxSIznmnlWaHQkMQk1MioiRYUdqE0uY2jGf7qEGJrQrVfkl/CE6NGJFr4MtKN+oD5ibmQtja552WmxoHDU1MTomIBkk3JP5vcx1sgu2jBk0GoGTbkLuk32FkroGmNhtWlpu/YlRs3ukOzjultOOrKBHRIK3a14KtdT5MLDI/och0Q/WSfgdNAG6HwrI9udClubEILR/KaISSjeYGQkMOk1MiokGQUuHtTbUwpEKuc2gmVMkSu6TfhKF6Sb/D2Hwde5odWFdjcmFU+7xTI7rP3DhoyGFySkQ0CF9VtGJdlRcTizjXdNCkf0hf0u+Q61CISoGP9ppbGCXae8wa0b2mxkFDD5NTIqIESanw1sYaRAyJAheb7g+WlK0AdAzVS/r7G52nY3WlG+Wt5v67EiIPMroTSkVNjYOGFianREQJWlPRirWVHpRw1HTQFCSU0QjAPqQv6XcYniPRErTh41KTR0+1fEjDw5ZSlFZMTomIEiClwhsbqmOjplyqdNCU9EKpEISwQAN6CxACKM6RWLYnF56QmW/VLgBhGNFSE2OgoYbJKRFRAlaWNeNrjpomTWy5UgkIXtLvMCZfR43fhk/3mdcFIrYMrwNGZBtbSlHaMDklIopT1JB4bV01DKWQz7mmg6YQhZTNABxmh2Ipdg1w2YElO/MQ1s2b6iC0Qki9qn1xBKLUY3JKRBSnf+1qxMYaLyYVs69pMiijBVCRIV+l35vxBTp2NzuwssK8EXoh8qBkAEZ0j2kx0NDC5JSIKA6+kI5Xvq6CQxPIcfAS9GApANJoACAAwbek7tx2BU0A/9ieb1pTfiE0QAgYkR3mBEBDDl8JiIji8PbmGuxpDGDSMI6aJoX0QUk/hODc3b5MKNSxtcGJLyvNHD0tgBHdBSX9psVAQweTUyKiAapoCeKtjTUoznHAYePLZzJIoxGAAQjO3e1LrkNBKuDtbfkwzBo91QqgpBdGdLc5AdCQwldXIqIBUErhb2sqUOeLYFwh2x0lg1IRSNkEFkId2MRCHZvqXFhl0uipEHYACnpkmynHp6GFySkR0QCsLGvBx7ubUFLshsYm8Ukhjcb2Qigm+weS64yNnv7vlgJEDXNiEKIg1lKKl/YpxZicEhEdgC+kY9HqcuhSoTiHo3zJoGBAGXUAbFwRaoBKinRsrndihUl9T4VWFLu0H9lpyvFp6GBySkR0AK+uq8K2Oj+mDGcRVLIoo6l9RSgWQg1UjiNWuf/3LQVoi6Y/oRftCyTokU1pPzYNLUxOiYj6sa7Sg8WbajA6zwkni6CSQkFCGrVg+6j4HVQcxc4mBz7YlWfK8YVWBCO6E9JgQ35KHb4qEBH1wRfS8fTKMgTCBkbls0F8siijGUoGOGqaAKcNyHcq/O/WfNT5099nV4gCKOmDEdma9mPT0MHklIioF0opLFpdjs21Phw8Ird9jXEarNioaU3sB8FFDBIxoVBHtc+Ov2/OR7qXuxdCA+BANLwGSpnU14qyHpNTIqJeLN/ThHe21GFcgQtOO18qk0UZTe1N9zl/N1GaAMbm61i6Jw8batPf6UDThkNGqyCje9N+bBoa+IpLRNRNeUsbnl65D0opjMjj5fxkUTAgjWrE5ppy1HQwhudIBHWBF9cXIqSnd1RfaDmACiMaXpfW49LQweSUiGg/gbCOPy7fi4rWIKaOyDU7nKwi9br2uaYcNR0sIYDJxTo21jmxeGt++o+vFcGIbGRhFKUEk1MionZSKjz9xT6sLm/FISNy2Ww/iZQKt881tbNCP0lcdoVit8SbW/OxvSG9/XeFVgRltEIPcfSUko+vEERE7V5fX4V/bK7FxCI33A5edk4WBUDqlQD7mibd2HwDrUEbnllTjEAkfR+mhNAgtFzo4VVQMpi249LQwOSUiAjARzsb8NevKlHotnMVqCRTsjW2VKlwczWoJBMCmDo8gg21Lry0oTCt1ftCGwFp1EPn3FNKMianRDTkfVXeisdXlEJKhXGFHNlLJgUdUi8HICEEi8tSwWUHxuTreGd7HpaXpW8+rxB2CDgRDX0GJUNpOy5lPyanRDSkbaz24pF/7YInpHN50iSLXc6vgJI+CMHislQakSuhhMCza4qwpzl9I//CNhJSr4Ie/jptx6Tsx+SUiIasjdVePLxsJ+r9EUwbyUb7yaZkE6Re3345n283qTa5OIo6vx2PfVGMprb0PN5COGKjp8HlUDKQlmNS9uOrBRENSWsrWvHg0h2o8YVx2Kg8JqZJplQQMroPgOLl/DTRBHDIiAg217vwp1XD0BZNz79pYRsNadQgEvw8Lcej7MfklIiGnE92N+JXy3ai0R/BYaPy2DIqyRR0GNG9UCrEy/lp5rABU4ZFsWJfDv6yuhgRI/XHFMIGIfKhh1ZA6rWpPyBlPSanRDRkSKnw+roq/Oaj3QiEDUxjYpp0ChIyWgolWyFEHqvzTZDrUCgpiuKDXbl4anUxoulIULURUEYrIoH3oVQaDkhZzW52AERE6eAL6XjmizK8u7UOBU47Jg1j8VOyKShIfR+k0RAbMeU8U9MUuhRUgYF3duRDAfjJcR647KnrMyWEgGYbAz2yEbbwV3C4T0jZsSj7MTkloqy3s96Px1eU4usqD0qK3OxjmgKxxLQ8dllXuADBtxezFbklAB3/2J6HYFRg/rdbUeBKYYKq5UJILyKBD2CzHwTNPj5lx6LsxlcPIspaEV3i7c01eGVtFRoDEUwbmQeXnaN5yaYgYyOmei0gnCyAspAit4RNKCzdk4eWkA23ntiC8QWpu+wubGMg9X0I+9+Cu/AaCI1zjil+Qql0rieRfF6vF0VFRfB4PCgsLDQ7HCKyiC21PrzwZTlWl7eg0GXHhCI3K/JTQEGHjJa2X8p3AUxMLSmkC+xtdmDKsCh+cnwrZo8Pp+xYSkUh9UrYc06EK/9iCI6iU5yYnBJRVqn1hvDG+mp8sL0B3nAUU4blItdpMzusrKRUAEa0FEp62ueYMgmxMkMCe5odyHEo/OfhPlw0w48cR2pSACXbIGUDHDlnwpk3F4LzjykOTE6JKCvU+8J4d2sd3t1Sh1pfGGMLnBiZ5+RoaQrE5pfWQRqVgIq0V+Uz+cgESgH1ARua22w4YmwYV87yYtbYcEqaKkjphZKtcOScAWfeeRCCHxJpYJicElHGUkphd2MAS3c04MOdjaj3hzEsx4FxhS62iEoBBQDSB0OvhJKtAGwQws12URkorANlrQ647QpnTAniwsN9OKhYT/pxOhJUu/vbcOXNg9DYJYMOjMkpEWWc1mAUq8tb8K9djdhQ7YU3pGNErgNj8l3QNCZKyaYAQAUg9VpIowmADiFyeBk/C7QENVT77BieY+C7U9pw9iFtOHh4NKmfN5QMQBoNsDkPhTP/AtjsE5O3c8pKTE6JyPKUUqj3R7Cx2os1Fa1YXd6KxkAYmhAYU+BCsdvOy/cpoCChpBfKaIA0WgDogHBBwMHR0iyiFNDQZkNjwIZCl8TR40I4c2oQx4wLIc+ZnBQhViRVBWErhMP9HThyTuYoKvXJEsnpk08+id/97neoqanBjBkz8Nhjj+HUU08d0O8yOSXKPkop1PrC2NvUhl0Nfqyt9KC0uQ2eYBQCAsNyHRiZ64DdxnmOyaYgAemHlB4oowlKBQEoJqVDgFJAa0hDnd8GIYAxeQZOmBjE0ePDOHxUBMNy5CD3r6BkE5TyQ7ONhyPnFNhdsyC0/CT9BZQtTE9OX3/9dVxxxRV48skncfLJJ+OZZ57Bf//3f2Pr1q046KCDDvj7TE6JMpdSCp6QjnpfGDXeEGp8YZQ2tWF7vR+NgQh8YR1SKrgdGoblOFDkdsDGy/ZJE7tcHwFUCEq1QUkfpPQCKgpAIjan1AWwkGXIiRhAY5sNnqANmqZQ7JaYOiyKmaPDmDIsioOKdIwr0OFI4J+GUjqU0QCFCIQ2DHbXUbA5psHmmMK+qATAAsnpt7/9bRxzzDF46qmnOm+bPn065s2bh4ULFx7w95mcElmLbkgEoxLBqIG2qIFAxIA/rMMX1uEN6fCEomj0R1DjDaPGF0IgbCAYNRDSJRQUNAjku2zId9lR4LLBrnF0NBGxF3YJKB2AAagoFHRARaBUBFBBKBmEQrR9GwVAALBDCAcTUuqky9iIqidkQ8QQEFDIcyjkOSUmFOqYVBzFmDwDw3IkinMMFLkkClwSBU6JHIdCX58nlTKgZCuU9AHCBqEVQLNNgOaYDJttFIRtOIRWBKEVsFfqEGPq2Y5EIli7di3uvvvuLrefffbZWLlyZa+/Ew6HEQ5/0zzY6/WmNMbe/O/aj/GvXdVpPy4NLfF8auz4iKm63PbNO4Lqts1gPpJKFft9qQCpRPv/AUMJGBIwFGDI2PdRKaBLQJcCUQOIGKLzNgBwaIDbruB2SLjtCh0R61GgNQq0+hOPM/MkelLUN/9X3X7u9UsCquPyrGhvAeWOfd9F6lYRokxkwKlF4dRiz3VfRENdwIatDS4IuOCyS7hsEk6bgtMm4bApODQFe/uXTVOwCQWbBmgilrDG/sUVt+9fto/YhwHsBMQeADbERu8FhNCgCdE+t1zr8q+1Y6ZJl3/B4ptbRM97qd3YAhfuPud7cNittaSzqclpY2MjDMPAmDFjutw+ZswY1NbW9vo7CxcuxAMPPJCO8Pq0ZFuVqcenLJZAftL3r8Tu6Z6Iqu7fq/3SmS7fi85EVnW7r7e9CQB2rf1Fpdugm10DnDbTp7dbXF+Pz3639/qpIo7HVdjAVaspORSArq2nlALChoBUPUfdFQBdAbrR9daOJFUIdI6waiL2CVhAxhLP9s9VonNPsZ0IiC5JKND9e7XfbbHvOWW6q7KWENoibSiyF5kdSheWeJXqXmWrlOqz8vaee+7B7bff3vmz1+tFSUlJSuPr7vnLL0/r8YiIiIiGClOT05EjR8Jms/UYJa2vr+8xmtrB5XLB5XKlIzwiIiIiSjNTKw2cTieOPfZYLFu2rMvty5Ytw0knnWRSVERERERkFtMv699+++244oorMHv2bJx44ol49tlnUV5ejhtuuMHs0IiIiIgozUxPTi+55BI0NTXhwQcfRE1NDWbOnIklS5Zg0qRJZodGRERERGlmep/TwWKfUyIiIqLswe7WRERERGQZTE6JiIiIyDKYnBIRERGRZTA5JSIiIiLLYHJKRERERJbB5JSIiIiILMP0PqeD1dEJy+v1mhwJEREREfWnoKAAQoh+t8n45NTn8wEASkpKTI6EiIiIiPozkL70Gd+EX0qJ6urqAWXiyeD1elFSUoKKigo2/c9QPIeZj+cw8/EcZjaev8xn1jkcEiOnmqZh4sSJaT9uYWEhn5AZjucw8/EcZj6ew8zG85f5rHgOWRBFRERERJbB5JSIiIiILIPJaZxcLhfuv/9+uFwus0OhBPEcZj6ew8zHc5jZeP4yn5XPYcYXRBERERFR9uDIKRERERFZBpNTIiIiIrIMJqdEREREZBlMTomIiIjIMpic9uLJJ5/ElClT4Ha7ceyxx2LFihX9br98+XIce+yxcLvdmDp1Kp5++uk0RUp9ieccvvXWW/i3f/s3jBo1CoWFhTjxxBPxz3/+M43RUm/ifR52+Pzzz2G32zFr1qzUBkj9ivf8hcNh/PznP8ekSZPgcrlw8MEH44UXXkhTtNSbeM/hyy+/jKOOOgq5ubkYN24crr76ajQ1NaUpWtrfp59+irlz52L8+PEQQuDtt98+4O9YKpdR1MVrr72mHA6Heu6559TWrVvVggULVF5entq3b1+v2+/du1fl5uaqBQsWqK1bt6rnnntOORwO9eabb6Y5cuoQ7zlcsGCBeuSRR9Tq1avVzp071T333KMcDof6+uuv0xw5dYj3HHZobW1VU6dOVWeffbY66qij0hMs9ZDI+Tv//PPVt7/9bbVs2TJVWlqqvvzyS/X555+nMWraX7zncMWKFUrTNPWnP/1J7d27V61YsULNmDFDzZs3L82Rk1JKLVmyRP385z9X//u//6sAqMWLF/e7vdVyGSan3Rx//PHqhhtu6HLbt771LXX33Xf3uv2dd96pvvWtb3W57frrr1cnnHBCymKk/sV7Dntz+OGHqwceeCDZodEAJXoOL7nkEvWLX/xC3X///UxOTRTv+Xv//fdVUVGRampqSkd4NADxnsPf/e53aurUqV1ue/zxx9XEiRNTFiMNzECSU6vlMrysv59IJIK1a9fi7LPP7nL72WefjZUrV/b6O1988UWP7c855xysWbMG0Wg0ZbFS7xI5h91JKeHz+TB8+PBUhEgHkOg5XLRoEfbs2YP7778/1SFSPxI5f++88w5mz56N3/72t5gwYQKmTZuGn/3sZwgGg+kImbpJ5ByedNJJqKysxJIlS6CUQl1dHd58803MmTMnHSHTIFktl7Gn/YgW1tjYCMMwMGbMmC63jxkzBrW1tb3+Tm1tba/b67qOxsZGjBs3LmXxUk+JnMPufv/73yMQCODiiy9ORYh0AImcw127duHuu+/GihUrYLfzZc1MiZy/vXv34rPPPoPb7cbixYvR2NiIG2+8Ec3NzZx3aoJEzuFJJ52El19+GZdccglCoRB0Xcf555+PP//5z+kImQbJarkMR057IYTo8rNSqsdtB9q+t9spfeI9hx1effVV/PKXv8Trr7+O0aNHpyo8GoCBnkPDMPDDH/4QDzzwAKZNm5au8OgA4nkOSikhhMDLL7+M448/Hueeey7+8Ic/4K9//StHT00UzzncunUrbrnlFtx3331Yu3YtPvjgA5SWluKGG25IR6iUBFbKZTjEsJ+RI0fCZrP1+GRYX1/f4xNFh7Fjx/a6vd1ux4gRI1IWK/UukXPY4fXXX8e1116Lv//97zjrrLNSGSb1I95z6PP5sGbNGqxbtw7z588HEEt2lFKw2+1YunQpzjzzzLTETok9B8eNG4cJEyagqKio87bp06dDKYXKykoceuihKY2ZukrkHC5cuBAnn3wy7rjjDgDAkUceiby8PJx66ql46KGHeBXR4qyWy3DkdD9OpxPHHnssli1b1uX2ZcuW4aSTTur1d0488cQe2y9duhSzZ8+Gw+FIWazUu0TOIRAbMb3qqqvwyiuvcI6UyeI9h4WFhdi0aRPWr1/f+XXDDTfgsMMOw/r16/Htb387XaETEnsOnnzyyaiurobf7++8befOndA0DRMnTkxpvNRTIuewra0NmtY1pbDZbAC+GYEj67JcLmNKGZaFdbTPeP7559XWrVvVrbfeqvLy8lRZWZlSSqm7775bXXHFFZ3bd7RfuO2229TWrVvV888/z1ZSJov3HL7yyivKbrerv/zlL6qmpqbzq7W11aw/YciL9xx2x2p9c8V7/nw+n5o4caK68MIL1ZYtW9Ty5cvVoYceqn70ox+Z9ScMefGew0WLFim73a6efPJJtWfPHvXZZ5+p2bNnq+OPP96sP2FI8/l8at26dWrdunUKgPrDH/6g1q1b19kKzOq5DJPTXvzlL39RkyZNUk6nUx1zzDFq+fLlnfddeeWV6vTTT++y/SeffKKOPvpo5XQ61eTJk9VTTz2V5oipu3jO4emnn64A9Pi68sor0x84dYr3ebg/Jqfmi/f8bdu2TZ111lkqJydHTZw4Ud1+++2qra0tzVHT/uI9h48//rg6/PDDVU5Ojho3bpy67LLLVGVlZZqjJqWU+vjjj/t9X7N6LiOU4ng7EREREVkD55wSERERkWUwOSUiIiIiy2BySkRERESWweSUiIiIiCyDySkRERERWQaTUyIiIiKyDCanRERERGQZTE6JiFLs8ccfhxACM2fO7HObvXv3Yv78+Zg2bRpycnKQm5uLGTNm4Be/+AWqqqp6bP+Pf/wDc+fOxZgxY+B0OjF8+HB897vfxcsvv4xoNJrKP4eIKKXYhJ+IKMVmzZqFDRs2AABWrVqFb3/7213uf/fdd3HppZdi5MiRmD9/Po4++mgIIbBp0ya88MIL0DQN69atAxBbp/yaa67BX//6V5x77rn44Q9/iJKSEng8Hnz88cd4/vnn8eCDD2LBggVp/zuJiJKBySkRUQqtWbMGxx13HObMmYP33nsP1113HZ599tnO+0tLS3HEEUdg2rRp+Pjjj1FUVNTl95VSWLx4Mb7//e8DAH7729/irrvuwgMPPID77ruvx/Fqa2uxe/dunHLKKan9w4iIUoTJKRFRCv3kJz/B008/jU2bNuH666/Hpk2bUFtbi9zcXADAzTffjCeeeAJffPEFTjjhhH73FY1GMXbsWIwePRpbt26FECIdfwIRUVpxzikRUYoEg0G8+uqrOO644zBz5kxcc8018Pl8+Pvf/965zdKlSzFmzJgDJqZAbBS2ubkZ3/ve95iYElHWYnJKRJQib775JjweD6699loAwCWXXIL8/Hw8//zznduUl5djypQpA9pfeXk5AAx4eyKiTMTklIgoRZ5//nnk5OTg0ksvBQDk5+fjoosuwooVK7Br1y6ToyMisiYmp0REKbB79258+umnmDNnDpRSaG1tRWtrKy688EIAwAsvvAAAOOigg1BaWjqgfR500EEAMODtiYgyEZNTIqIUeOGFF6CUwptvvolhw4Z1fs2ZMwcA8OKLL8IwDJxzzjmoq6vDqlWrDrjP2bNnY/jw4fi///s/sJaViLIVk1MioiQzDAMvvvgiDj74YHz88cc9vn7605+ipqYG77//Pm677Tbk5eXhxhtvhMfj6bGvjlZSAOBwOHDXXXdh+/bt+NWvftXrsevr6/H555+n9O8jIkoltpIiIkqyd999F3PnzsUjjzyCO++8s8f9jY2NmDhxIv7jP/4DixcvxrvvvotLLrkEo0eP7mzCDwBbt27tHIHtrQn/nDlzujTh//TTT/Hss8/igQceYBN+IspYTE6JiJLsggsuwJIlS1BZWYlRo0b1us0PfvADvPnmm6isrMSYMWOwd+9e/P73v8fSpUtRUVEBTdMwZcoU/Pu//ztuvvlmTJ48ucvvv/POO3j22WexevVqtLS0oKCgALNmzcIll1yCq6++Gk6nMw1/KRFR8jE5JSIiIiLL4JxTIiIiIrIMJqdEREREZBlMTomIiIjIMpicEhEREZFlMDklIiIiIstgckpERERElsHklIiIiIgsg8kpEREREVkGk1MiIiIisgwmp0RERERkGUxOiYiIiMgymJwSERERkWX8PyOYg9KWhvPqAAAAAElFTkSuQmCC","text/plain":["<Figure size 800x400 with 1 Axes>"]},"metadata":{},"output_type":"display_data"}],"source":["import numpy as np\n","import matplotlib.pyplot as plt\n","from scipy.stats import beta, binom\n","import scipy.stats as st\n","from matplotlib.lines import Line2D  \n","import seaborn as sns\n","\n","# 定义正确率的取值范围\n","acc_values = np.linspace(0, 1, 1000)\n","\n","# 定义先验分布 Beta(70, 30)\n","alpha_prior = 70\n","beta_prior = 30\n","prior_distribution = beta(alpha_prior, beta_prior).pdf(acc_values)\n","\n","# 定义似然分布 Bin(50, ACC)\n","n_trials = 50\n","y_observed = 30\n","likelihood = binom.pmf(y_observed, n_trials, acc_values)\n","\n","# 绘图\n","fig, ax = plt.subplots(figsize=(8, 4))\n","\n","ax.fill_between(acc_values, prior_distribution, color=\"#f0e442\", alpha=0.6, label='prior')\n","likelihood_scaled = likelihood / np.max(likelihood) * np.max(prior_distribution)  # 缩放\n","ax.fill_between(acc_values, likelihood_scaled, color=\"#0071b2\", alpha=0.6, label='(scaled) likelihood')\n","\n","# 图例和标题\n","ax.set_xlabel('ACC', fontsize=12)\n","ax.set_ylabel('Density', fontsize=12)\n","ax.legend(loc='upper right', fontsize=10)\n","\n","# 创建自定义图例\n","custom_lines = [Line2D([], [], color=\"#f0e442\"),\n","                Line2D([], [], color=\"#0071b2\")]\n","        \n","# 移除图的上、右边框线\n","sns.despine()"]},{"cell_type":"markdown","id":"35bd50d1","metadata":{"id":"290B46A55E204F4D80D18AFF8E7E41A6","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["**思考🧐**  \n","\n","哪一张图正确反映了正确率ACC的后验模型？  \n","\n","![Image Name](https://cdn.kesci.com/upload/skg2tifhs.png?imageView2/0/w/960/h/960)\n"]},{"cell_type":"markdown","id":"2f08e6fe","metadata":{"id":"56553C9245874C0A8190FE0399E1B104","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["我们可以画出后验分布图："]},{"cell_type":"code","execution_count":10,"id":"fe9ed856","metadata":{"collapsed":false,"id":"1400BA605E734E23A004CBD9EF7F44B4","jupyter":{"outputs_hidden":false},"notebookId":"66ea43e4f99628c505715aea","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/1400BA605E734E23A004CBD9EF7F44B4/sjznude5kl.png\">"],"text/plain":["<Figure size 1000x400 with 1 Axes>"]},"metadata":{},"output_type":"display_data"}],"source":["# 导入统计建模工具包 scipy.stats 为 st\n","import scipy.stats as st \n","\n","# 设置 x 轴范围 [0,1]\n","x = np.linspace(0,1,10000)\n","# 设置 Beta 分布参数\n","a,b = 70,30\n","# 形成先验分布 \n","prior = st.beta.pdf(x,a,b)/np.sum(st.beta.pdf(x,a,b))\n","\n","# 形成似然\n","k = 30     # k 代表正确率为1的次数\n","n = 50     # n 代表总次数\n","likelihood = st.binom.pmf(k,n,x)\n","\n","# 计算后验\n","unnorm_posterior = prior * likelihood                  # 计算分子\n","posterior = unnorm_posterior/np.sum(unnorm_posterior)  # 结合分母进行计算\n","likelihood = likelihood /np.sum(likelihood)            # 为了方便可视化，对似然进行类似后验的归一化操作 \n","\n","# 绘图\n","plt.plot(x,posterior, color=\"#009e74\", alpha=0.5, label=\"posterior\")\n","plt.plot(x,likelihood, color=\"#0071b2\", alpha=0.5, label=\"likelihood\")\n","plt.plot(x,prior, color=\"#f0e442\", alpha=0.5, label=\"prior\")\n","plt.legend()\n","plt.xlabel(\"ACC\")\n","plt.fill_between(x, prior, color=\"#f0e442\", alpha=0.5)\n","plt.fill_between(x, likelihood, color=\"#0071b2\", alpha=0.5)\n","plt.fill_between(x, posterior, color=\"#009e74\", alpha=0.5)\n","sns.despine()"]},{"cell_type":"code","execution_count":11,"id":"2a13fa5c","metadata":{"collapsed":false,"id":"C7E84E0382BC48049034760576331C1D","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"outputs":[],"source":["# ------"]},{"cell_type":"markdown","id":"1fc45905","metadata":{"id":"318579D7884A4A04829C28501422F255","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"subslide"},"tags":[]},"source":["**正式计算**  \n","\n","$$  \n","P(\\text{ACC} | y = 30) = \\frac{P(\\text{ACC})L(\\text{ACC}|y = 30)}{P(y = 30)}.  \n","$$  \n","\n","和之前一样，分母$P(y = 30)$是一个常数，在计算中可以将其忽略  \n","\n","\n","$$  \n","P(\\text{ACC} | y = 30) \\propto \\left[ \\frac{\\Gamma(100)}{\\Gamma(70)\\Gamma(30)} \\cdot \\binom{50}{30} \\right] \\cdot \\text{ACC}^{69} (1 - \\text{ACC})^{29} \\cdot \\text{ACC}^{30} (1 - \\text{ACC})^{20}  \n","$$  \n","\n","<center>[  ] 中的也是可以忽略的常数项<center>  \n","\n","\n"]},{"cell_type":"markdown","id":"8e602273","metadata":{"id":"623D6CD69E7A4ED395C5D58A1D73C901","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"subslide"},"tags":[]},"source":["整理一下可知，后验分布可以表示为：  \n","$$  \n","P(\\text{ACC} | y = 30) \\propto \\text{ACC}^{99} (1 - \\text{ACC})^{49}  \n","$$  \n","\n","\n","根据这个公式，我们发现 $P(\\pi | y=30)$ 和 $Beta(100,50)$ 有着相似的形状  \n","\n","$$  \n","\\text{Beta}(100, 50) = \\frac{\\Gamma(150)}{\\Gamma(100)\\Gamma(50)} \\pi^{99} (1-\\pi)^{49}  \n","$$  \n","实际上，在这里后验分布也确实是Beta分布：  \n","\n","$$  \n","\\pi | (Y = 30) \\sim \\text{Beta}(100,50)  \n","$$"]},{"cell_type":"markdown","id":"9385d064","metadata":{"id":"1B4E10488C7A42A38732870E9D9FECA1","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["**对后验模型进行总结**  \n","\n","在结合先验和似然之后，我们对正确率ACC的认识发生了更新。  \n","\n","需要注意的是：后验模型仍然是一个$Beta$分布，和先验模型是同一个分布族。  \n","\n","$$  \n","ACC | (Y = 30) \\sim \\text{Beta}(100,50)  \n","$$  \n","\n","$$  \n","f(\\text{ACC} | y = 30) = \\frac{\\Gamma(200)}{\\Gamma(100)\\Gamma(100)} \\text{ACC}^{99} (1-\\text{ACC})^{99}  \n","$$  \n"]},{"cell_type":"markdown","id":"5f69400a","metadata":{"id":"E41DF1343C074EBDA6349466B3B8F93D","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["所以我们也可以对二者进行对比，下表进行了这一总结，可以发现，在新的数据产生之后：  \n","\n","例如：  \n","\n","- 对正确率的期望值从0.70降低为0.67；  \n","\n","- 模型的标准差从0.0657减少为0.0471；  \n","\t\n","|    |prior  |posterior  \n","|----|-----|----|  \n","|$\\alpha$   |70  |100 |  \n","|$\\beta$   |30  |50 |  \n","|mean   |0.70  |0.67 |  \n","|mode  |0.694  |0.663 |  \n","|var   |0.00432  |0.00222 |  \n","|sd   |0.0657  |0.0471 |  "]},{"cell_type":"markdown","id":"33effbab","metadata":{"id":"09A819A33DFB45D69C1082C91CD1EA97","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["**思考🧐**  \n","\n","- 先验分布的形态是否决定了后验分布的形态？(还记得之前关于正态分布和Beta分布的例子吗)  \n","- 先验分布和后验分布是否必然一致？"]},{"cell_type":"markdown","id":"8c0f3fdb","metadata":{"id":"568799DE04984CDDAE10AE5F6DA0E829","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["## Part 4: Beta-Binomial model and simulation"]},{"cell_type":"markdown","id":"061926c1","metadata":{"id":"FB88FC13435B467EAACD048F6437AE6E","jupyter":{"outputs_hidden":false},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["**Beta-Binomial model 的一般数学形式**  \n","\n","在前一节中，我们对随机点运动任务的正确率ACC建立了基本、更加完整的Beta-Binomial模型。  \n","\n","- 我们现在假设了一个特定的Beta(45,55)先验分布和特定的正确判断次数（正确次数： Y = 30，总试次： n = 50）。  \n","- 这只是Beta-Binomial模型的一个特例。 实际上 Beta 模型可以适用于任何参数范围在[0,1]的场景。  \n","- 例如，硬币为正面的倾向，可重复研究的比率或是疫苗的有效率。  \n","\n","$$  \n","\\begin{split}  \n","Y | ACC & \\sim \\text{Bin}(n, ACC) \\\\  \n","ACC & \\sim \\text{Beta}(\\alpha, \\beta). \\\\  \n","\\end{split}  \n","$$  \n"]},{"cell_type":"markdown","id":"92fcc51e","metadata":{"id":"EEC81CF8269242B2A72E2F404E6E8A80","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["无论是哪种情况，在观察到 n 次时间中有 $Y = y$次目标事件后，ACC的后验分布可以用Beta模型来描述，反映了先验（通过$\\alpha$和$\\beta$）和数据（通过$y$和$n$）的影响：  \n","\n","$$  \n","\\begin{equation}  \n","ACC | (Y = y) \\sim \\text{Beta}(\\alpha + y, \\beta + n - y)  .  \n","\\end{equation}  \n","$$  \n","\n","\n","需要注意的是，后验与先验是相同的概率模型，只是参数不同  \n","\n","- 在这个例子中，Beta($\\alpha$, $\\beta$)模型是对应数据模型的 $\\text{Bin}(n, ACC)$ 的**共轭先验 (conjugate prior)**。  \n","- 如果 $f(ACC)$ 是 $L(ACC|y)$ 的共轭先验，后验 $f(ACC|y) \\propto f(ACC)L(ACC|y)$ 与先验来自相同的模型族。  \n","- 我们将在第五次课详细介绍共轭先验的相关知识。  \n"]},{"cell_type":"markdown","id":"6804ab19","metadata":{"id":"F83423A889E54387AD966E49987E915F","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["**代码实现**  \n","\n","让我们模拟正确率ACC的后验模型。  \n","\n","- 首先，我们从Beta(45,55)先验中模拟10,000个ACC值  \n","- 然后, 使用从每个ACC值中模拟Bin(50,ACC)的潜在正确判断次数Y：  \n","\n","模拟结果：10,000对ACC和y值的模拟数据  \n"]},{"cell_type":"code","execution_count":12,"id":"eb223f2d","metadata":{"collapsed":false,"id":"DA95A13370E24F5381A45C6CE43FEB6D","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>ACC</th>\n","      <th>y</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>0</th>\n","      <td>0.479492</td>\n","      <td>28</td>\n","    </tr>\n","    <tr>\n","      <th>1</th>\n","      <td>0.385205</td>\n","      <td>16</td>\n","    </tr>\n","    <tr>\n","      <th>2</th>\n","      <td>0.412269</td>\n","      <td>19</td>\n","    </tr>\n","    <tr>\n","      <th>3</th>\n","      <td>0.458158</td>\n","      <td>20</td>\n","    </tr>\n","    <tr>\n","      <th>4</th>\n","      <td>0.397845</td>\n","      <td>24</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["        ACC   y\n","0  0.479492  28\n","1  0.385205  16\n","2  0.412269  19\n","3  0.458158  20\n","4  0.397845  24"]},"execution_count":12,"metadata":{},"output_type":"execute_result"}],"source":["# 导入数据加载和处理包：pandas\n","import pandas as pd\n","# 导入数字和向量处理包：numpy\n","import numpy as np\n","# 导入基本绘图工具：matplotlib\n","import matplotlib.pyplot as plt\n","\n","# 设置随机种子，以便后续可以重复结果\n","np.random.seed(84735)\n","\n","# 模拟 10000 次数据\n","king_sim = pd.DataFrame({'ACC': np.random.beta(45, 55, size=10000)})  # 从Beta(45,55)先验中模拟10,000个ACC值\n","king_sim['y'] = np.random.binomial(n=50, p=king_sim['ACC'])       # 从每个ACC值中模拟Bin(50,ACC)的潜在正确判断次数Y\n","\n","# 显示部分数据\n","king_sim.head()"]},{"cell_type":"markdown","id":"dec5f579","metadata":{"id":"1F0BCD78B26444A59F7DA3EB9F55051B","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["通过散点图观察以上数据的关系:  \n","- 黑色点表示正确次数$Y\\neq30$的部分。  \n","- 蓝色点表示正确次数$Y=30$的部分。"]},{"cell_type":"code","execution_count":13,"id":"e32ca1bf","metadata":{"collapsed":false,"id":"C8E64540304A4B08AD062610297933DD","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/C8E64540304A4B08AD062610297933DD/sjznudppnz.png\">"],"text/plain":["<Figure size 1000x400 with 1 Axes>"]},"metadata":{},"output_type":"display_data"}],"source":["# 绘制散点图：正确次数 (Y!=30)部分，用黑色表示\n","plt.scatter(king_sim['ACC'][king_sim['y']!=30], \n","            king_sim['y'][king_sim['y']!=30], \n","            c='black', s = 3,\n","            label='FALSE')\n","# 绘制散点图：正确次数 (Y=30)部分，用蓝色表示\n","plt.scatter(king_sim['ACC'][king_sim['y']==30],\n","            king_sim['y'][king_sim['y']==30],\n","            c='b', s = 20, \n","            label='TRUE')\n","\n","# 显示图片\n","plt.legend(title = \"y==30\")\n","plt.xlabel('ACC')\n","plt.ylabel('Y')\n","plt.show()"]},{"cell_type":"markdown","id":"9583b959","metadata":{"id":"C622D4345F044151A3B09AA25DD66E2C","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["当我们仅关注与我们的Y = 30正确次数匹配的对时，剩下的ACC值的很好地逼近了后验模型Beta(75,75)："]},{"cell_type":"code","execution_count":14,"id":"0b53c017","metadata":{"collapsed":false,"id":"3F77076BCD944A41B8A9B779E89A13DD","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/3F77076BCD944A41B8A9B779E89A13DD/sjznuezofg.png\">"],"text/plain":["<Figure size 500x500 with 1 Axes>"]},"metadata":{},"output_type":"display_data"}],"source":["# 导入绘图工具 seaborn 为 sns\n","import seaborn as sns\n","\n","# 从模拟数据中筛选出 y 值为 30 的样本，生成对应的后验分布。\n","king_posterior = king_sim[king_sim['y'] == 30]\n","\n","# 绘制分布图：概率密度+柱状图\n","sns.displot(king_posterior['ACC'], kde=True)\n","plt.xlabel('ACC')\n","plt.ylabel('Density')\n","plt.show()"]},{"cell_type":"markdown","id":"22878240","metadata":{"id":"5D76A021D70F4C4F91AE2E9866E3F744","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["同时绘制出先验、似然和后验。"]},{"cell_type":"code","execution_count":15,"id":"c13a06cf","metadata":{"collapsed":false,"id":"ECA4D86EE7A94D3D9AAD899B8BAAF74A","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/ECA4D86EE7A94D3D9AAD899B8BAAF74A/sjznuetdtz.png\">"],"text/plain":["<Figure size 1000x400 with 1 Axes>"]},"metadata":{},"output_type":"display_data"}],"source":["# 导入统计建模工具包 scipy.stats 为 st\n","import scipy.stats as st \n","\n","# 设置 x 轴范围 [0,1]\n","x = np.linspace(0,1,10000)\n","# 设置 Beta 分布参数\n","a,b = 45,55\n","\n","# 形成先验分布 \n","prior = st.beta.pdf(x,a,b)/np.sum(st.beta.pdf(x,a,b))\n","\n","# 形成似然\n","k = 30                  # k 代表正确率为1的次数\n","n = 50                  # n 代表总次数\n","likelihood = st.binom.pmf(k,n,x)\n","\n","# 计算后验\n","unnorm_posterior = prior * likelihood                  # 计算分子\n","posterior = unnorm_posterior/np.sum(unnorm_posterior)  # 结合分母进行计算\n","likelihood = likelihood /np.sum(likelihood)            # 为了方便可视化，对似然进行类似后验的归一化操作 \n","\n","# 绘图\n","plt.plot(x,posterior, color=\"#009e74\", alpha=0.5, label=\"posterior\")\n","plt.plot(x,likelihood, color=\"#0071b2\", alpha=0.5, label=\"likelihood\")\n","plt.plot(x,prior, color=\"#f0e442\", alpha=0.5, label=\"prior\")\n","plt.legend()\n","plt.xlabel(\"ACC\")\n","plt.fill_between(x, prior, color=\"#f0e442\", alpha=0.5)\n","plt.fill_between(x, likelihood, color=\"#0071b2\", alpha=0.5)\n","plt.fill_between(x, posterior, color=\"#009e74\", alpha=0.5)\n","plt.xlim([0,1])\n","plt.show()"]},{"cell_type":"code","execution_count":16,"id":"b0329e9b","metadata":{"collapsed":false,"id":"A8CB201DB2B64C03B9AA573CE02589D4","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"outputs":[],"source":["# ------"]},{"cell_type":"markdown","id":"8c053b98","metadata":{"id":"7CDE7DF37857491EBDB0D0DEB9733E48","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["我们还可以使用模拟样本近似后验特征，例如随机点运动任务正确率的均值和标准差。  \n","\n","- 结果与上面计算的理论值非常相似，$E(ACC∣Y = 30) = 0.5024$ 和 $SD(ACC∣Y = 30) = 0.0394$：  \n","- 在解释这些模拟结果时，“近似”是一个关键词。由于我们10,000次模拟中只有219次与观测到的Y = 30数据匹配，通过将原始模拟次数从10,000增加到50,000，可以改善这个近似值："]},{"cell_type":"code","execution_count":17,"id":"51b0e5e5","metadata":{"collapsed":false,"id":"013D5BA7004049D1AE206B29EBB1BF82","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"outputs":[{"name":"stdout","output_type":"stream","text":["近似值： 均值， 0.5024772340319912 。标准差， 0.0394983348061003\n"]}],"source":["print(\n","  \"近似值：\", \n","  \"均值，\",\n","  king_posterior['ACC'].mean(), \n","  \"。标准差，\",\n","  king_posterior['ACC'].std()\n",")"]},{"cell_type":"code","execution_count":18,"id":"ede114d8","metadata":{"collapsed":false,"id":"21458B0299E94816926D4F3463E79BBD","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"outputs":[{"name":"stdout","output_type":"stream","text":["10,000次模拟中, 219次与观测到的Y = 30数据匹配\n","50,000次模拟中, 1045次与观测到的Y = 30数据匹配\n"]}],"source":["print(f\"10,000次模拟中, {king_posterior.shape[0]}次与观测到的Y = 30数据匹配\")\n","\n","# 模拟新的数据\n","size = 50000 # 不同于之前的 10000\n","king_sim2 = pd.DataFrame({'ACC': np.random.beta(45, 55, size=size)})\n","king_sim2['y'] = np.random.binomial(n=50, p=king_sim2['ACC'])\n","king_posterior2 = king_sim2[king_sim2['y'] == 30]\n","print(f\"50,000次模拟中, {king_posterior2.shape[0]}次与观测到的Y = 30数据匹配\")"]},{"cell_type":"markdown","id":"cf76b245","metadata":{"id":"0B05C83508AE4E6AB5EA1DE8E717C746","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["### <center> 🎈总结 </center>  \n","\n","在第3章中，我们对一个真实的现象(随机点运动任务中特定条件下正确判断的能力)，构建了Beta-Binomial模型，并通过数据进行了一次信念更新：  \n","\n","$$  \n","\\begin{split}  \n","Y | ACC &\\sim \\text{Bin}(n, ACC) \\\\  \n","ACC &\\sim \\text{Beta}(\\alpha, \\beta)  \n","\\end{split}  \n","\\Rightarrow ACC | (Y = y) \\sim \\text{Beta}(\\alpha + y, \\beta + n - y).  \n","$$  \n"]},{"cell_type":"markdown","id":"cd6bc4d7","metadata":{"id":"7C826D66F04148929A51F37A89FD1802","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["这个模型反映了贝叶斯分析的四个通用要素:  \n","\n","1. 先验模型 Beta先验模型可以通过调节参数来反映不同的ACC值在[0,1]范围内的相对先验可能性。  \n","\n","$$  \n","f(ACC) = \\frac{\\Gamma(\\alpha + \\beta)}{\\Gamma(\\alpha)\\Gamma(\\beta)}ACC^{\\alpha - 1}(1-ACC)^{\\beta - 1}.  \n","$$  \n","\n","2. 数据模型 为了学习ACC,我们收集数据$Y$,即$n$次独立试验中成功的次数,其中每次试验正确的概率为ACC。$Y$对ACC的依赖关系由二项分布$Bin(n, ACC)$描述。  \n","\n","3. 似然函数 在观测到数据$Y=y$,其中$y\\in{0,1,...,n}$后,ACC的似然函数通过将$y$代入二项式概率质量函数而获得,它提供了一种机制来比较不同ACC与数据的兼容性:  \n","\n","$$  \n","L(ACC|y) = \\binom{n}{y}ACC^y(1-ACC)^{n-y} \\text{ for } ACC \\in [0,1].  \n","$$  \n","\n","4. 后验模型 通过贝叶斯规则,共轭的Beta先验和二项式数据模型结合产生ACC的Beta后验模型。更新的Beta后验参数($\\alpha$ + $y$, $\\beta$ + $n$ - $y$)反映了先验的影响(通过$\\alpha$和$\\beta$)和观测数据的影响(通过$y$和$n$)。  \n","\n","$$  \n","f(ACC|y) \\propto f(ACC)L(ACC|y) \\propto ACC^{(\\alpha+y)-1}(1-ACC^{(\\beta+n-y)-1}.  \n","$$  \n"]},{"cell_type":"markdown","id":"1f98f2cb","metadata":{"id":"3582752C620C463085298360D4B9FA64","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["**🥋练习**  \n","\n","- 假设心理学研究中的可重复研究出现的概率为 0.4（这里概率用$\\pi$表示）， 请自行选择一个合适的 Beta 分布，模拟10000个$\\pi$值。  \n","- 请根据$\\pi$模拟相应的数据 $Y$ (可重复研究的数量)，假设总的可重复研究数量 $n = 100$。  \n","- 绘制先验分布和后验分布的图像。"]},{"cell_type":"markdown","id":"ff7d2274","metadata":{"id":"ADCB20A2829748E3A7837AC0880312D6","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["### 💐Bonus 1：Beta(70,30)这个先验是怎么选取的?"]},{"cell_type":"markdown","id":"67ef5886","metadata":{"id":"DC41A41FE9114B3EA1C5816D96470E6F","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["**Beta分布的集中趋势量数**  \n","\n","在回答这个问题之前，我们先来了解一下**Beta分布的集中趋势量数**  \n","\n","需要注意的是，虽然我们可能比较熟悉计算正态分布的均值、众数和标准差；但 Beta 分布的相关统计量的计算有所不同：  \n","\n","![Image Name](https://cdn.kesci.com/upload/s0ywp66o5v.png?imageView2/0/w/500/h/500)  \n",">source:https://en.wikipedia.org/wiki/Probability_density_function  \n","\n","**1. 平均数(mean)**  \n","\n","- ACC的平均取值  \n","$$  \n","E(ACC)  = \\frac{\\alpha}{\\alpha + \\beta}  \n","$$  \n","$$  \n","E(ACC) = \\int ACC \\cdot f(ACC)dACC.  \n","$$  \n","\n","**2. 众数(mode)**  \n","\n","- ACC最可能的取值。即，在ACC下，$f(ACC)$能取到最大值。  \n","\n","$$  \n","\\text{Mode}(ACC)  = \\frac{\\alpha - 1}{\\alpha + \\beta - 2} \\;\\;\\; \\text{ when } \\; \\alpha, \\beta > 1.  \n","$$  \n","$$  \n","\\text{Mode}(ACC) = \\text{argmax}_{ACC}f(ACC).  \n","$$  \n","**3. 方差(variance)**  \n","\n","- $ACC$取值的变异性(variability)  \n","$$  \n","\\text{Var}(ACC) = E((ACC - E(ACC))^2) = \\int (ACC - E(ACC))^2 \\cdot f(ACC) dACC.  \n","$$  \n","\n","$$  \n","\\text{Var}(ACC) = \\frac{\\alpha \\beta}{(\\alpha + \\beta)^2(\\alpha + \\beta + 1)} .  \n","$$  \n","\n","**4. 标准差(standard deviation)**  \n","$$  \n","\\text{SD}(ACC) = \\sqrt{\\text{Var}(ACC)} .  \n","$$  \n","\n","-------  \n","-------------"]},{"cell_type":"markdown","id":"41b52569","metadata":{"id":"C11F18C0B1D9498BB8874B066C3B8022","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["**调整Beta先验**  \n","\n","我们已经知道，在先前的推测中，在 5% 一致性条件下，被试的平均正确率约为 70%。根据这一点，我们可以计算出 Beta 分布的$\\alpha$和 $\\beta$参数。  \n","$$  \n","E(ACC) = \\frac{\\alpha}{\\alpha + \\beta} = 0.70  \n","$$  \n","重新整理过后：  \n","$$  \n","\\alpha = 0.70(\\alpha + \\beta)  \n","$$  \n","\n","$$  \n","\\alpha \\approx \\frac{7}{3} \\beta  \n","$$  \n","\n","一个合适的分布中，$\\alpha$和$\\beta$需要满足的条件如上，我们可以选择 *Beta(70,30), Beta(7,3), Beta(14,6)*，我们可以通过以下代码来画出这些分布的形状"]},{"cell_type":"code","execution_count":19,"id":"95ee1e85","metadata":{"collapsed":false,"id":"97C43C38CAE84331963D8D7EEBB87186","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/97C43C38CAE84331963D8D7EEBB87186/sjznueox6r.png\">"],"text/plain":["<Figure size 1300x500 with 1 Axes>"]},"metadata":{},"output_type":"display_data"}],"source":["import numpy as np\n","from scipy.stats import beta\n","import matplotlib.pyplot as plt\n","\n","fig, ax = plt.subplots(figsize=(13, 5))\n","pz.Beta(70, 30).plot_pdf(pointinterval=True,ax=ax)\n","pz.Beta(7, 3).plot_pdf(pointinterval=True,ax=ax)\n","pz.Beta(14, 6).plot_pdf(pointinterval=True,ax=ax)\n","plt.tight_layout()\n","plt.show()"]},{"cell_type":"markdown","metadata":{"id":"DACC8775EFC94756BA0F8594A6F91169","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["如图所示，我们选择$ACC \\sim \\text{Beta}(70,30)$作为先验模型  \n","\n","带入公式，可以计算出$先验f(ACC)$的pdf、平均数、众数、方差和标准差：  \n","\n","- pdf:  \n","$$  \n","f(ACC) = \\frac{\\Gamma(100)}{\\Gamma(70)\\Gamma(30)}ACC^{69}(1-ACC)^{29} \\;\\; \\text{ for } ACC \\in [0,1].  \n","$$  \n","\n","- 平均数:  \n","$$  \n","E(ACC) = \\frac{70}{70 + 30} = 0.70  \n","$$  \n","\n","- 众数:  \n","$$  \n","\\text{Mode}(ACC) = \\frac{70 - 1}{70 + 30 - 2} = 0.6939  \n","$$  \n","\n","- 方差:  \n","$$  \n","\\text{Var}(ACC) = \\frac{70 \\cdot 30}{(70 + 30)^2(70 + 30 + 1)} = 0.0021  \n","$$  \n","\n","- 标准差：  \n","$$  \n","\\text{SD}(ACC) = \\sqrt{0.0021} = 0.0458  \n","$$"]},{"cell_type":"markdown","metadata":{"id":"C3E7DF2F45224164ACD047E25642951B","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["### 💐Bonus2 Beta-Binomial 模型生成的后验仍是Beta分布"]},{"cell_type":"markdown","metadata":{"id":"2A9DDD921BDB46DD857FC54578EBBA6C","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["**公式推导**  \n","\n","对于 $ACC$ 的后验模型为 $Beta(α + y, β + n − y)$ 的推导过程 (其中，正确次数为y，总试次为n)。  \n","\n","1. 首先写出先验的公式  \n","\n","$$  \n","\\begin{equation}  \n","f(ACC) = \\frac{\\Gamma(\\alpha + \\beta)}{\\Gamma(\\alpha)\\Gamma(\\beta)}ACC^{\\alpha - 1}(1-ACC)^{\\beta - 1}  \n","\\;\\; \\text{ and } \\;\\;  \n","L(ACC|y) = \\left(\\!\\!\\begin{array}{c} n \\\\ y\\end{array}\\!\\!\\right) ACC^{y} (1-ACC)^{n-y}  .  \n","\\end{equation}  \n","$$  \n","\n","2. 结合先验和似然函数 (暂时忽略分母)，后验概率密度函数可以由贝叶斯定理得到：  \n","\n","$$  \n","\\begin{split}  \n","f(ACC | y)  \n","& \\propto f(ACC)L(ACC|y) \\\\  \n","& = \\frac{\\Gamma(\\alpha + \\beta)}{\\Gamma(\\alpha)\\Gamma(\\beta)}ACC^{\\alpha - 1}(1-ACC)^{\\beta - 1}  \\cdot \\left(\\!\\begin{array}{c} n \\\\ y \\end{array}\\!\\right) ACC^{y} (1-ACC)^{n-y} \\\\  \n","& \\propto ACC^{(\\alpha + y) - 1} (1-ACC)^{(\\beta + n - y) - 1}  .\\\\  \n","\\end{split}  \n","$$  \n"]},{"cell_type":"markdown","metadata":{"id":"6BB243523A8849CB9C643C46ADDF5805","jupyter":{},"notebookId":"66ea43e4f99628c505715aea","runtime":{"execution_status":null,"is_visible":false,"status":"default"},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":["\n","3. 最后，我们加上归一化因子(分母部分)：  \n","\n","$$  \n","f(ACC|y) = \\frac{\\Gamma(\\alpha+\\beta+n)}{\\Gamma(\\alpha+y)\\Gamma(\\beta+n-y)}ACC^{(\\alpha + y) - 1} (1-ACC)^{(\\beta + n - y) - 1}.  \n","$$  \n","\n","我们知道$Beta(a,b)$的概率密度函数可以写为：  \n","$$  \n","Beta(α, β) = \\frac{\\Gamma(\\alpha + \\beta)}{\\Gamma(\\alpha)\\Gamma(\\beta)} ACC^{\\alpha-1} (1-ACC)^{\\beta-1} \\;\\;  for\\; ACC \\in [0,1] \\  \n","$$  \n","\n","因此，  \n","\n","$$  \n","f(ACC|y) = Beta(α + y, β + n − y)  \n","$$"]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.12.4"}},"nbformat":4,"nbformat_minor":5}
